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INTRODUCTION 


This  final  report  covers  the  period  18  February  1966 
through  1  March  1967,  the  period  of  performance  under  contract 
AF  33(657) -15919.  The  purpose  is  to  accumulate,  in  one  report, 
a  concise  discussion  of  the  research,  the  operations,  and  the 
capabilities  of  the  Seismic  Data  Laboratory  (SDL)  as  they  relate 
to  the  VELA  UNIFORM  Program. 

The  Seismic  Data  Laboratory  was  originally  established 
under  Contract  AF  33 (657) -7427  and  was  operated  as  the  Data 
Analysis  and  Technique  Development  Center  from  October  1961 
through  16  August  1963.  Basically,  this  contract  period  was 
devoted  to  the  design  and  establishment  of  a  permanent  facility 
which  would  serve  not  only  as  a  repository  for  all  the  data  accu¬ 
mulated  for  use  under  the  VELA  UNIFORM  Program,  but  also  as  an 
organization  to  analyze  the  data,  conduct  basic  research  in  seismic 
phenomena,  and  to  perform  data  services  for  the  entire  VELA  UNIFORM 
community. 

These  goals  were  essentially  completed  in  February  1963 
when  the  SDL  moved  into  its  present  facility  at  300  N.  Washington 
Street,  Alexandria,  Virginia.  By  this  time,  the  organization 
structure  was  established,  personnel  with  the  qualifications 
needed  to  perform  the  work  were  on  the  job,  routine  operating 
procedures  had  been  established,  the  electronic  data  processing 
equipment  configuration  was  complete,  and  a  research  plan  had  been 
established  to  attack  the  problem  of  the  detection  and  identifica¬ 
tion  of  earthquakes  and  nuclear  explosions. 

Thus  the  objectives  of  this  earlier  peiiod  had  been  attained, 
and  the  Seismic  Data  Laboratory  began  its  operations  as  an  integral 
part  of  the  VELA  UNIFORM  Program  under  contract  AF  33(657) -12447. 
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The  work  performed  and  the  accomplishments  attained  under 
this  contract  were  reported  in  the  final  report  issued  on  15  April 

1966. 

Essentially,  the  work  performed  under  the  present  contract 
was  a  continuation  of  the  tasks  outlined  under  prior  SDL  contract 
periods.  The  actual  work  statement  is  as  follows: 

1.  Tasks 

a.  Perform  data  analyses  and  seismological  research  as 
directed  by  the  AFTAC  project  officer.  The  following  task  assign¬ 
ments  are  included: 

(1)  Design,  develop,  and  apply  new  analytical  techniques 

and  procedures  for  processing  seismological  data. 

(2)  Analyze  and  evaluate  data  from  explosions  and 

selected  earthquakes. 

(3)  Compile  and  report  basic  data  from  explosions  and 
earthquakes,  which  may  include  the  variatior  of  amplitude  with  dis¬ 
tance  as  a  function  of  frequency,  travel  times  of  body  and  surface 
waves,  travel-time  anomalies,  variation  of  amplitude  with  yield  oc 
explosion  environment,  wave-form  characteristics,  and  other  infor¬ 
mation  as  specified. 

(4)  Evaluate  existing  identification  criteria  and 
develop  and  apply  new  identi  .ication  criteria  for  distinguishing 
between  seismic  waves  from  natural  and  artificial  sources. 

(5)  Develop  data  processing  techniques  that  wij.1  im¬ 
prove  the  detection  and  identification  capability  of  station 

networks  and  associated  analysis  centers. 

(6)  Process  seismological  data  from  the  Large  Aperture 

Seismic  Array  as  directed  by  the  AFTAC  project  officer. 
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b.  Maintain  data  services  and  plant  facilities  as 

follows: 

(1)  Maintain  a  central  data  file,  to  include  a  system 
for  cataloging  and  classifying  incoming  seismic  data,  that  include 

means  for  rapid  access  to  significant  data. 

(2)  Maintain  a  file  of  all  digitized  seismograms. 

(3)  Continually  evaluate  the  equipment,  facilities, 
and  staff  of  the  Seismic  Data  Laboratory  (SDL);  make  recommenda¬ 
tions  for  modifications  as  required,  and  effect  such  modifications 

as  approved  by  the  Government. 

(4)  Obtain  and  prepare  digital  and  analog  computer 

programs  for  cataloging,  processing,  and  analyzing  the  data. 

(5)  Furnish  copies  of  seismological  recordings  or 
processed  seismological  data  to  other  '£LA  UNIFORM  participants, 
provide  facilities  for  visiting  scientists  to  review  data  within 
the  SDL,  provide  computational  services  to  other  VELA  participants, 
conduct  mathematical  analyses,  and  other  activities  as  approved  by 
the  AFTAC  project  officer. 

(6)  Compress  the  data  on  analog  magnetic  tapes  so  that 
the  original  tapes  can  be  degaussed  and  returned  to  the  field  for 
further  use  as  directed  hy  the  AFTAC  project  officer. 

These  tasks  essentially  cover  the  work  which  has  been  per¬ 
formed  during  this  period  and  the  summary  discussion  which  follows 
is  based  on  all  the  detailed  reports  and  memoranda  which  have  been 
issued  under  the  contract. 

Section  II  discusses  the  general  SDL  facilities,  \ncluding 
analog  and  digital  machine  data  processing.  Section  III  is  con¬ 
cerned  with  the  storage  and  retrieval  of  film  and  tape,  tape  com¬ 
pression,  and  the  data  which  have  been  processed  in-house.  Section 
IV  covers  the  types  of  service  functions  which  SDL  offers  to  other 
VELA  UNIFORM  participants.  Section  V  is  a  descriptive  discussion 
of  the  statistical  summary  analyses  completed  on  earthquake  and  ex 
plosion  events,  and  Section  VI  gives  a  concise  summary  of  the  re¬ 
search  studiec  which  have  been  undertaken. 


*  its 
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11 *  FACILITIES 

A.  Office  Lay-Out 

The  Seismic  Data  Laboratory  is  housed  in  the  second  and 
third  floors  of  a  seven-  story  building  at  300  North  Washington  Street, 
Alexandria,  Virginia.  A  floor  plan  of  these  two  floors  is  shown  in 
Figure  1.  In  general  the  offices  are  on  the  third  floor  while  the 
second  floor  is  restricted  to  computer  activities  and  the  Central 
Data  Files.  Owing  to  security  requirements,  normal  access  to  the 
laboratory  is  entirely  through  the  third  floor  lobby  door  that  is 
constantly  attended  by  a  receptionist  during  the  daytime  shift. 

Other  exits  are  provided  !nto  the  stairweJis  but  can  be  used  only 
in  case  of  emergency  by  breaking  a  security  seal.  Also,  as  needed, 
freight  can  be  entered  either  through  an  outside  freight  door  in 
the  maintenance  room  or  the  second  floor  elevator  door  that  is  nor¬ 
mally  kept  locked. 

In  addition  to  the  regular  SDL  staff,  office  space  is 
provided  for  the  AFTAC  Project  Officer,  for  other  AFTAC  personnel, 
and  for  SDL  visitors.  This  provision  fulfills  one  of  the  functions 
of  SDL.  Visitors,  either  VELA  UNIFORM  participants  or  approved 
foreign  scientists,  are  encouraged  to  use  the  facilities  provided 
tc  study  the  data  stored  in  SDL. 

The  present  building  includes  many  features  required 
specifically  for  the  computing  facilities.  These  features  include 
oxtra  electric  power,  separate  high-capacity  air  conditioners,  and 
a  raised  floor  to  accommodate  wiring  and  to  serve  as  an  air  con¬ 
ditioning  plenum  under  the  computer  room  and  adjacent  maintenance 
area.  The  regular  concrete  flooring  under  this  area  has,  in  ad¬ 
dition,  extra  reinforcing  steel  for  special  loading. 
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igure  1.  Seismic  Data  Laboratory  Floor  Plan 
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The  Central  Data  Files  provides  apace  for  250  three- 
foot  shelf  sections,  each  of  which  holds  102  magnetic  tapes.  The 
total  storage  capacity  is  approximately  25,550  tapes.  A  large 
portion  of  data  over  one  year  old  is  being  compressed,  thus  in¬ 
creasing  greatly  the  actual  amount  of  data  that  can  eventually 
be  stored  here.  The  compressed  data  tapes  are  stored  at  another 
near-by  facility  in  order  that  the  most  current  data  will  be 

readily  available  for  processing. 

Approximately  300  square  feet  is  also  devoted  to 

storing  analog  film  at  SDL.  New  techniques  for  storing  digital 
tapes  have  enabled  us  to  store  all  our  digital  tapes  in  the  com¬ 
puting  area. 

B.  Personnel  and  Management 

The  SDL  is  organized  into  three  sections  in  accordance 
with  Figure  2.  The  Data  Analysis  Section  consists  of  trained  geo¬ 
physicists  who  make  routine  examination  of  seismic  records?  the 
Research  Section  consists  of  scientists  on  the  Ph.D.  level  and  is 
responsible  for  the  original  research  conducted;  and  the  Computa¬ 
tional  Services  Section  provides  the  electronic  data  processing 
function.  As  a  secondary  function,  the  Data  Analysis  and  Compu¬ 
tational  Service-,  sections  perform  service  for  other  VELA  UNIFORM 
participants. 

The  execution  of  a  research  program  is  ideally  per¬ 


formed  on  a  project  basis  rather  than  an  organization  basis.  The 
project  leader  is  one  of  the  senior  research  scientists  supported 
by  one  or  more  computational  assistants,  analysts  or  programmers. 
Since  a  research  worker  cannot  keep  the  same  number  of  support 
people  busy  all  the  time,  pooling  the  various  specialties  makes 
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for  more  efficient  use  of  personnel.  The  use  of  service  pools  j.s 
also  an  efficient  way  of  providing  necessary  support  for  non-SDL 
requests  without  untenable  interruption  of  research  programs.  This 

applies  to  both  analysts  and  programmers. 

The  SDL  has  made  use  of  consultants  to  design  experiments 
to  follow  their  progross,  and  to  help  interpret  the  results;  the 
actual  work  being  done  by  geophysicists  or  senior  analysts  assigned 
full  time  to  the  projects  and  supervision  being  furnished  by  appro¬ 


priate  SDL  personnel . 

C.  Electronic  Data  Processing  Equipment 

The  equipment  configuration  is  indicated  schematically 
in  Figures  3  and  4.  Thi3  equipment  may  be  roughly  classified  into 
functional  groups  as  follows;  (1)  magnetic  tape  playback,  including 
filters,  spectrum  analyzers,  etc.,  (2)  analog  computer,  (3)  time 
code  generator  and  the  anaiog-to-digital  converter,  (4)  electronics 
for  compressing  data  on  magnetic  tape,  and  (5)  digital  computer  and 
its  supporting  equipment. 

The  various  modes  of  data  flow  through  the  digital  com¬ 
puter  are  shown  schematically  in  Figure  3.  The  most  convenient 
aspects  of  this  system  are:  (1)  The  small  CDC  160-A  that  performs 
offline  operations  freeing  the  1604  computer  for  the  bigger  jobs, 
and  (2)  the  plotter  that  has  been  used  almost  universally  to  present 
the  output  in  a  form  readily  understood  by  visual  inspection.  Plot¬ 
ting  will  now  be  shared  by  the  recently  completed  D/A  conversion 
system,  including  the  analog  plotter. 


In  order  to  machine  process  LASA  data  as  set  forth  in 
the  work  statement,  the  equipment  configuration  was  modified  to 
include  the  following  additions: 

■  One  CDC  604  magnetic  tape  drive 

-  One  CDC  818  disk  file 

-  One  CDC  162-3  controller 

-  A  digital-to-analog  sub-system 

The  disk  file,  tape  drive,  and  the  controller  were 
supplied  and  installed  by  the  Control  Data  Corporation.  Because 
of  our  specialized  requirements,  the  individual  parts  of  the  D/A 
sub-system  were  purchased,  assembled  and  made  operational  by  SDL 
engineering  personnel. 

The  analog  sections  of  the  equipment  are  characterized 
by  high  flexibility.  Any  element  can  be  connected  through  patch- 
panels  to  any  other  element.  The  various  elements  in  the  system 
are  shown  in  their  relative  floor  arrangement  in  Figure  4;  the  con¬ 
necting  points  for  the  trunk  cables  between  the  patch  panels  are 
as  shown. 

The  following  operations  can  now  be  performed,  many 
simultaneously,  with  the  analog  processing  system: 

1.  Transcribe  the  magnetic  recordings  from  one  tape  to 
another,  either  to  make  a  duplicate  tape  or  to  produce  a  conposite 
tape  consisting  of  selected  portions  of  many  original  tapes.  Whereas 
the  original  tape  is  recorded  at  0.3  ips,  the  transcription  can  be 
performed  at  speeds  up  to  60  ips,  that  is,  200  times  faster; 

2.  Make  paper  records,  or  film  playouts,  from  magnetic 

recordings; 
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3.  Make  power  spectra  analysis  and  plot  the  results 
automatically; 

4.  Digitise  analog  tapes  for  further  analysis  on  the 
digital  computer.  The  digital  tapes  are  in  an  IBM  compatible 
format; 

5.  Generate  and  read  VELA  Time  Code; 

6.  Filter  data  using  Krohn-Hite  filters  or  special 
analog  computer  circuits; 

7.  Make  develocorder  films  from  analog  tape; 

8.  Read  and  write  analog  tapes  compatible  with  the 
British  24-channel  system; 

9.  Compress  data  from  field  tapes  on  new  tapes  in 
the  ratio  of  10:1,  thus  freeing  nine  tapes  out  of  ten  for  re-use 
in  the  Field  Measurement  program; 

10.  Transcribe  data  from  standard  IRIG  Jj-inch  magnetic 
tapes  to  the  one-inch  tapes  normally  used  at  SDL;  and 

Perform  any  general  analog  computer  operations. 
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Ill .  DATA  STORAGE  AND  RETRIEVAL 

A.  Magnetic  Tape  and  Film 

Seismic  data  comes  to  the  SDL  from  the  Long  Range  Seismic 
Measurements  (LRSM)  field  teams,  U.S.  and  foreign  observatories,  a  d 
from  the  Large  Aperture  Seismic  Array  (LASA)  located  at  Billings, 
Montana.  The  FM  analog  magnetic  tapes  and  operational  logs  are  stored 
in  the  Central  Data  Files  and  the  digital  tapes  from  LASA  are  stored 
in  the  digital  computer  area.  For  working  convenience,  both  16  mm 
and  35  mm  film,  together  with  operational  logs,  are  stored  adjacent 
to  the  Data  Analysis  Section. 

When  this  data  comes  in,  it  is  checked  against  a  packing 
list  and  a  receiving  record  is  prepared.  This  receiving  record  con¬ 
tains  the  recording  station's  name  and  the  date  of  recording.  The 
data  itself  is  stored  by  station  and  month.  When  data  is  needed  for 
processing,  the  receiving  record  is  checked  to  see  if  the  data  has 
been  received  and,  if  so,  the  cataloger  goes  to  the  appropriate 
shelving  unit  where  the  data  is  sotied  and  checks  it  out  of  the  files 
for  processing. 

B.  Data  Tape  Compression  Program 

Early  in  1963  it  became  clearly  evident  that  the  storage 
of  the  analog  magnetic  tape  would  soon  be  a  storage  problem  if  all 
data  recorded  in  the  field  was  to  be  stored  in  the  SDL  facility. 
Because  of  this,  and  to  partially  solve  the  high  cost  of  magnetic 
tapes  used  in  the  field  program,  it  was  decided  to  compress  this 
original  data  by  re-recording  it  on  magnetic  tape  at  10  tijnes  its 
present  density;  that  is  to  say,  the  data  on  10  original  tapes 
would  be  recorded  on  one  new  tape. 
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Certain  seismic  data,  particularly  that  of  known  nuclear 
explosions,  is  maintained  in  its  original  form.  Originally,  data 
six  months  old  was  compressed,  but  once  the  initial  backlog  of 
data  was  eliminated,  only  data  at  least  one  year  old  is  compressed. 

The  basic  system  requires  the  use  of  three  tape  transports 
that  may  be  manipulated  by  a  single  operator.  While  the  data  is 
being  transferred  from  the  tape  on  one  transport  to  the  tape  on 
another  transport,  a  third  transport  is  being  loaded  and  prepared 
for  transcription. 

As  the  original  data  is  being  compressed,  a  log  is  pre¬ 
pared  which  contains  the  recording  station  and  data  day,  the  data 
compressed,  quality  control  check-off,  and  the  date  the  tape  is 
erased.  The  date  that  the  data  is  compressed  is  noted  in  the 
proper  receiving  record,  discussed  above.  The  box  containing  the 
compressed  data  is  marked  with  the  station  name  and  data  recording 
day.  The  compressed  tapes  are  then  stored  at  another  near-by 
facility  where  they  are  readily  available. 

Compression  operations  started  in  September  1963,  an'? 
the  sub-system  became  fully  operational  during  February  1964. 

Present  production  is  being  maintained  at  the  level  needed  to  meet 
the  requirements  of  the  field  measurement  program  (LRSM  teams  and 
U.S.  Observatories)  and  SDL. 

C.  Composite  Tapes 

The  data  library  contains  analog  composite  tapes  which 
were  either  prepared  by  the  SDL  or  obtained  from  other  contractors 
for  use  in  the  VELA  UNIFORM  program.  A  separate  area  is  maintained 


-10- 


within  the  Central  Data  Files  for  storage  of  these  composite  topes 
and  log-books  are  kept  by  both  the  Central  Data  Files  and  Analog 
Processing  supervisor.  These  log-books  contain  the  event  name,  and 
the  recording  stations  which  are  a  part  of  the  composite  tape.  A 
composite  tape  for  a  given  event  contains  each  station  calibration, 

10  minutes  of  data  recorded  at  each  station  before  the  P  arrival, 
and  10-20  minutes  of  data  recorded  after  the  P  arrival. 

D.  Digitized  Data 

The  digital  library  contains  many  thousands  of  digital 
seismograms,  most  of  which  were  recorded  on  analog  tape  and  digi¬ 
tized  at  SDL.  Several  hundred  events  were  recorded  directly  on 
digital  tape  uexng  data  acquisition  systems  at  Tonto  Forest  Se  .s  no- 
logical  Observatory  and  LASA.  In  each  case  data  have  been  de¬ 
multiplexed  into  a  library  format  which  can  be  input  to  digital 
programs  already  in  use. 

Digitized  data  can  be  punched  on  cards  or  written  on 
tape  in  an  IBM  compatible  format.  The  data  is  written  in  the  binary 
mode  at  200  or  556  bits  per  inch  (as  specified  by  the  requestor) . 

The  tapes  contain  a  variable  number  of  files  of  36-bit  word  binary 
integers.  Each  file  contains  either  a  vertical  or  horizontal  com¬ 
ponent  from  a  single  station  for  a  specified  period  of  time  (up  to 
13  minutes) . 

The  binary  integers  range  between  plus  and  minus  2048. 

The  data  is  not  necessarily  calibrated,  but  instructions  are  avail¬ 
able  for  relating  digital  counts  to  earth  motion  at  leps. 

•"he  sampling  rate  can  be  varied,  but  is  normally  20  samples 
per  second.  The  data  is  prefiltered  to  minimize  aliasing.  The  time 
of  the  first  digital  point  can  be  determined  to  +  .1  second. 
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A  separate  list  is  furnished  with  each  tape  giving  the 
order  of  the  seismograms  written  on  each  tape,  sampling  rate,  number 
of  files  per  seismogram,  number  of  data  points  per  file,  time  of 
first  point,  component  identification,  event  name,  recording  station, 
and  filter  setting  used  at  time  of  digitization. 

After  each  seismogram  is  written  in  library  format,  it  is 
plotted  on  an  incremental  x-y  plotter  and  sent  through  a  QC  proce¬ 
dure.  If  the  QC  is  satisfactory,  magnitude,  time  of  first  point, 
time  of  first  motion,  signal- to-noise  ratio,  channel  quality,  event 
information  and  instrument  or  configuration  data  is  punched  on 
cards  together  with  the  location  on  magnetic  tape  of  each  seismogram. 
The  digitized  data  can  then  be  searched  on  any  of  the  above  param¬ 
eters  to  determine  its  availability  or  suitability  for  a  given 
analysis  assignment. 

E.  Digital  Programs 

Many  of  the  scientific  and  systems  programs  which  have 
either  been  developed  by  the  SDL  or  obtained  from  other  sources  have 
been  written-up  and  are  on  punched  cards  or  tape,  or  both. 

F.  U.S.  Coast  &  Geodetic  Survey  Epicenters 

BCD  tapes  containing  epicenters  listed  on  the  Preliminary 
Determination  of  Epicenters  ( PDE)  cards  (January  1960  to  date)  are 
n  file.  Geographic  and  seismic  regions  have  been  computed  for 
each  epicenter. 

'Pie  cards  contain: 

Orig.n  Time  Magnitude 

Latitude  Source  of  Magnitude 

Longitude  Geographic  Region  Number 

Depth  (KM)  Seismic  Reqion  Numbei 

Comments 
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G.  Earthquake  Bulletin  Data 


Earthquake  bulletin  data  from  the  LRSM  teams  (February 
1962  to  date)  and  the  VELA  observatories  (February  1963  to  date) 
hive  been  written  on  magnetic  tape.  Data  includes  the  PDE  cards, 
plus  all  phase  information  for  those  stations  which  recorded  the 
earthquake.  Recorded  phase  arrivals  not  associated  with  an  epi¬ 
center  are  also  on  the  tape. 


H.  Shot  Report  Data 


Data  from  more  than  75  shot  events  have  been  punched 
on  cards.  The  data  on  the  punched  cards  includes 

Shot  name  First  motion  (if  picked) 


Location 


Period 


Magnitude 


Amplitude 


Medium 


Instrument 


Recording  station 


Component 


Phase  ID 


Azimuth 


Arrival  Time 


Distance 


IV.  SERVICE  FUNCTIONS 


A.  Computational  Services  for  Other  VELA  Participants 

Manpower  and  oo^uter  time  spent  in  providing  services 
to  other  VELA  participants  has  increased  steadily  since  the  inception 
of  the  SDL  project.  During  this  contract  period,  286  requests  were 
received  from  50  organizations  and  the  VELA  Seism*  logical  Center. 
Table  1  lists  the  organizations  which  received  SDL  services  during 
the  period  of  this  report. 

The  normal  servic  s  3  .VLA  participants  include  copies  of 
16  and  35  mm  film,  copies  >f  digital  programs,  digitized  data,  copies 
of  analog  tape  composites,  and  the  use  of  the  1604  digital  computer 
for  checking  out  new  programs  or  running  production  programs. 
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ORGANIZATIONS  RECEIVING  SDL  DATA  SgRyXCg 


Advanced  Research  Projects  Agency 
AFCRL 

Bolt,  Beranek,  &  Newman,  Inc. 

Boston  College 

Brown  Engineering,  Inc. 

California  Institute  of  Technology 
Canadian  Department  of  Mines 
Canadian  Pacific  Oil  &  Gas  Company 
Carnegie  Institute  of  Washington 
Continental  Oil  Company 
Control  Data  Corporation 

Dominion  Observatory 
Department  of  Commerce 

Earth  Sciences,  Teledyne 
Engineering  Physics 

General  Atronics  Corporation 
Geotech,  Teledyne 
Graduate  Research  Center  of 
the  Southwest 

IBM  Corporation 

Institute  for  Defense  Analysis 
International  Seismological  Research 
Centre 

Larnont  Geophysical  Observatory 
Lawrence  Radiation  Laboratory 
Lincoln  Laboratory 

Massachusetts  Institute  of  Technology 
Mitre  Corporation 


Roland  F.  Beers,  Inc. 

Rowe  Air  Force  Base 

San  Calixto,  Bolivia  Observatory 
Stanford  Research  Institute 
Stanford  University 
St.  Louis  University 

Texas  Instruments,  Incorporated 
Triangle  Institute 

Underwater  Systems,  Inc, 

United  Kingdom  Atomic  Energy 
Authority 

University  of  Bergen 
University  of  Frankfort , Germany 
University  of  Michigan 
University  of  Minnesota 
University  of  Tasmania 
U.  S.  Coast  and  Geodetic  Survey, 
ESSA 

Vitro  Corporation 

Washington  Science  Center  (USCGS) 

Xavier  Univeysity 


Oregon  State  University 


Pan  American  Petroleum 
Pennsylvania  State  University 
Fhilco  Corporation 
Princeton  University 


Table  1 


As  the  seismic  community  has  become  familiar  with  the 
activities  of  the  SDL  and  aware  of  the  data  which  is  available, 
more  frequent  requests  have  been  received  for  the  following  services: 

-  Digitizing  data  from  specified  earthquakes  or  explosions, 

-  Running  SDL  production  programs,  such  as  power  spectral 
density  and  the  linear  array  programs,  on  specified  data. 

-  Writing  digitized  data  in  specified  formats  for  use  on 
other  computers. 

-  Computing  lists  of  earthquakes  which  meet  specified  criteria 
of  depth,  magnitude,  and  location. 

-  Making  composite  analog  tapes  of  special  events. 

-  Providing  space  and  assistance  for  visiting  scientists  to 
study  data  and  exchange  information  with  the  SDL  personnel. 

Data  requests  from  VELA  participants  are  screened  care¬ 
fully  so  as  not  to  interfere  with  normal  SDL  analysis  assignments, 
and  the  services  performed  are  approved  by  the  VSC. 

B.  Automated  Bulletin  Process 

Under  VSC' a  direction,  an  Automated  Bulletin  Process 
(ABP)  computer  program  was  developed  jointly  by  SDL  and  the  Geotech 
Division,  Teledyne  Industries,  Inc.  The  program  associates  seismic 
phase  arrivals  at  a  number  of  stations,  with  epicenters  reported  by 
the  U.S.  Coast  and  Geodetic  Survey,  and  identifies  up  to  23  phases. 
The  ABP  is  now  being  used  to  accelerate  the  preparation  of  the 
monthly  earthquake  bulletins  of  the  five  VELA  UNIFORM  seismological 
observatories  published  by  the  Geotech  Division,  Teledyne  Industries, 
Inc. 

In  terms  of  time  and  manpower  requirements,  most  of  the 
work  of  preparing  the  earthquake  bulletins  involves: 
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1.  Measurement  of  arrival  times,  amplitude,  and  period, 
preliminary  phase  identification  based  on  signal  character,  and 
card-punching  these  data  for  each  station.  Analysts  may  record 
from  300  to  2,000  events  per  month  at  each  station,  depending 
upon  the  station,  time  of  year,  and  seismic  activity  for  a  given 
month. 

2.  Association  of  phase  arrivals  with  hypocenters  re¬ 
ported  in  the  U.S.  Coast  and  Geodetic  Survey  (USC&GS)  Preliminary 
Determination  of  Epicenters  ( PDE)  cards,  identification  of  phases, 
merging  of  associated  and  unassociated  events,  computation  of 
azimuth,  distance,  magnitude,  and  ground  displacement,  and  finally, 
combination  of  all  data  from  the  five  observatories  into  a  form 
suitable  for  publication. 

The  ABP  was  designed  to  automate  the  latter  task. 

The  ABP  has  been  written  in  three  parts  which  can  be 
combined  into  one  program  but  are  usually  run  independently. 

V.  DESCRIPTIVE  ANALYSES  OF  EARTHQUAKES  AND  EXPLOSIONS 

A  long  range  seismic  measurements  (LRSM)  program  was  estab¬ 
lished  under  VELA  UNIFORM  to  record  and  analyze  short-period  and 
long-period  data  from  a  planned  series  of  U.S.  underground  nuclear 
teats.  A  standard  format  was  established  for  compiling  atatiscical 
data  from  these  nuclear  explosions  and  certain  earthquakes  having 
special  characteristics.  After  an  event  was  analyzed,  reports 
were  issued  summarizing  these  statistics  to  be  used  by  VELA  UNIFORM 
par  icipanta  in  studying  and  developing  methods  for  distinguishing 
between  explosive  and  earthquake  sources. 

The  data  from  these  various  events  were  rt corded  by  operating 
mobile  field  teams  anu  the  observatories  operated  at  Wichita  Moun¬ 
tain,  Oklahoma  (WMSo) ,  Uinta  Basin,  Utah  (UBSO) ,  Blue  Mountain, 
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Oregon  (BMfiO) ,  Cumberland  Plateau,  Tennessee  (CPSO) ,  and  Tonto 
Forest,  Arizona  (TFSO) ,  and  from  several  experimental  or  temporary 
stations  operated  in  connection  with  other  research  programs. 

Instrumentation  at  each  of  the  mobile  stetiors  consists  of 
three-component  short-period  Benioff  and  three-component  Sprengnether 
long-period  seismographs.  Data  are  recorded  on  35  millimeter  film 
and  one-inch  14  channel  magnetic  tape.  All  of  these  stations  are 
equipped  to  record  WWV  continuously  in  order  to  provide  accurate 
time  control.  Calibration  is  accomplished  once  each  day  and  just 
prior  to  each  shot  at  operating  settings.  Specific  details  of  the 
instrumentation  and  operating  procedures  for  these  statiors  are 
given  in  Field  Manual,  Long  Range  Seismic  Measurement  Program, 
Technical  Report  No.  63-17,  which  can  be  obtained  from  the  Geotech 
Division  of  Toledyne  Industries,  Inc.,  Dallas,  Texas.  All  the 
observatories  have  both  long-period  and  short-period,  three- 
component  instrumentation  in  addition  to  their  other  specialized 
facilities. 

Station  site  information,  presented  in  these  reports,  includes 
the  station  name  and  code;  the  geographic  coordinates,  distances 
and  azimuths  involved;  the  station  elevations;  and  the  type  of  in¬ 
struments  in  use  at  each  location. 

For  each  event  analyzed,  a  status  report  is  included,  giving 
the  names  of  stations  and  indicating  which  instruments  were  opera¬ 
tional  and  which  recorded  signals,  together  with  an  operations  map 
showing  the  location  of  each  station  recording  the  event  and  the 
type  of  signals  recorded. 


An  explanation  of  the  procedure  for  amplitude  measurement® 
used  in  these  reports  is  illustrated  in  Figure  5  .  The  unified 
magnitude  (m)  computations  for  distances  less  than  16°  are  based 
on  AFTAC/VSC  extensions  of  Gutenberg's  Tables*.  For  these  purposes, 
points  from  10°  to  16°  were  read  from  a  curve  in  the  Gutenberg- 
Richter  paper  ani  an  inverse  cube  relationship  was  used  to  extrap¬ 
olate  from  two  to  ten  degrees.  A  table  of  the  distance  factors  (B) 
is  shown  in  Table  2. 

A  standard  hypocenter  location  program  for  a  digital  computer 
is  used  to  determine  the  event  locations,  using  data  from  all  sta¬ 
tions  analyzed.  Best-fit  values  of  latitude,  longitude,  depth  of 
focus,  and  t' me  of  origin  were  determined  statistically  by  a  least 
squares  technique.  This  utilized  a  Jeffreys-Bullen  travel-time 
curve  as  modified  by  Herrin  in  1961  on  the  basis  of  Pacific  surface- 
focus  recordings.  Precision  of  the  computation  is  limited  primarily 
by  the  accuracy  of  arrival  times,  the  validity  of  the  standard 
travel-time  curve,  and  by  local  velocity  deviations.  Since  the 
method  is  based  on  P  wave  arrivals,  this  particular  program  does 
not  make  use  of  later  phases  such  as  pP  and  S  in  the  determination 
of  depth  or  location. 

The  measurements  made  of  the  principal  phases  from  these 
events  are  summarized.  These  include  the  Pn  and  P  arrival  times, 
the  maximum  amplitudes  (A/T)  of  Pn,  P  and  Pg  motion  as  seen  on  the 
short-period  vertical  instruments,  and  the  maximum  amplitudes  (A/T) 
of  the  Lg  phase  as  measured  on  the  short-period  horizontal  tan¬ 
gential  component.  Long-period  Rayleigh  wave  motions  are  also 
tabulated  in  (A/T)  form. 


♦Gutenberg,  B.  and  Richter,  C.F. ,  Magnitude  and  Energy  of  Earth 
quakes,  Ann.  Geofis.,  9  (1956),  pp.  1-15. 


Maximum 


Bottom  of  line  I 


Detail  Showing  Allowance 
For  Line  Width 


Pick  time  of  Pn  at  beginning  of  "a"  half  cycle. 

Pick  amplitude  of  Pn  as  maximum  "d/^"  within  2  or  3  cycles  of  "e". 
Pick  amplitudes  of  Pg  and  Lg  at  maximum  of  corresponding  motion. 


Figure  5.  Amplitude  Measurement*!  Procedures  for  Statistical 

Summary  Reports 
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Unified  Magnitude:  m  -  log1Q  ^A//T^'  4  B 

where  A  ■  zero  to  peak  ground  motion  in  millimicrons 

-  (mm)  (1000) 

K 

T  *  signal  period  in  seconds 
B  *  distance  factor  (see  Table  below) 
mm  »  record  amplitude  in  millimeters  zero  to 
peak 

K  •  magnification  in  thousands  at  signal 
frequency 

Table  of  Distance  Factors  (B)  for  Zero  Depth 


Dist 

Dist 

Dist 

Dist 

Idea).. 

_ B 

l.deg)_ 

B 

(deq) 

_ B 

(deg) 

0° 

— 

27° 

3.5 

54° 

3.8 

80° 

3.7 

1 

2 

3 

4 

2.2 

2.7 

3.1 

28 

29 

30 

31 

3.6 

3.6 

3.6 

3.7 

55 

56 

57 

58 

3.8 

3.8 

3.8 

3.8 

81 

82 

83 

84 

3.8 

3.9 
4.0 
4.0 

5 

3.4 

32 

3.7 

59 

3.8 

85 

4.0 

6 

7 

8 

9 

3.6 

3.8 

4.0 

4.2 

33 

34 

35 

36 

3.7 

3.7 

3.7 

3.6 

60 

61 

62 

63 

3.8 

3.9 

4.0 

3.9 

86 

87 

88 

89 

3.9 

4.0 

4.1 

4.0 

10 

4.3 

37 

3.5 

64 

4.0 

90 

4.0 

11 

12 

13 

14 

4.2 

4.1 

4.0 

3.6 

38 

39 

40 

41 

3.5 

3.4 

3.4 

3.5 

65 

66 

67 

68 

4.0 

4.0 

4.0 

4.0 

91 

92 

93 

94 

4.1 

4.1 

4.2 
4.1 

15 

3.3 

42 

3.5 

69 

4.0 

95 

4.2 

16 

17 

18 

19 

2.9 

2.9 

2.9 

3.0 

43 

44 

45 

46 

3.5 

3.5 

3.7 

3.8 

70 

71 

72 

73 

3.9 

3.9 

3.9 

3.9 

96 

97 

98 

99 

4.3 

4.4 

4.5 
4.5 

20 

3.0 

47 

3.9 

74 

3.8 

100 

4.4 

21 

22 

3.1 

3.2 

46 

49 

3.9 

3.8 

75 

76 

3.8 

3  .9 

101 

102 

4.3 

4.4 

23 

24 

3.3 

3.3 

50 

51 

3.7 

3.7 

f  U 

77 

78 

3.9 

3.9 

103 

104 

4.5 

4.6 

» 

25 

3.5 

52 

3.7 

79 

3.8 

105 

4  7 

26 

3.4 

53 

3.7 

Table  2.  Unified  Magnitudes  From  PR  or  P  Waves 
for  Statistical  Summary  Reports 


Other  data  given  in  theae  reports  include  the  travel-time 
residuals  from  the  Pn  and  P  phases  and  the  amplitudes  of  Pn  and  P, 
Pg  and  Lg,  (which  also  showed  lines  proportional  to  the  inverse 
cube  of  the  distance  visually  fitted  to  the  observed  points) , 
Rayleigh  wave  amplitudes,  and  illustrative  seismograms  showing 
the  signals  recorded  at  a  number  of  locations. 

Table  3  gives  the  events  analysed  by  SDL  together  with  basic 
data  for  each  event. 


Duff 

16  December  1965 

Tuff 

5.14  + 

0.51 

Charcoal 

10  September  1965 

Tuff 

5.16  + 

0.29 

Chartreuse 

6  May  1966 

Rhyolite 

5.22  1 

0.62 

Dumont 

19  May  1966 

Tuff 

'>.48  + 

0.56 

buryca 

14  April  1966 

Rhyolite 

5.17  + 

0.50 

Half  Beak 

30  Jun«»  1966 

Rhyolite 

6.02  + 

0.60 

Palanquin 

14  April  1965 

Rhyolite 

4.33  + 

0.36 

Pin  Stripe 

25  April  1966 

Tuff 

4.51  2. 

0.56 

Redhct 

5  March  1966 

Tuff 

4.11  + 

0.30 

Rex 

24  February  1966 

Tuff 

4.80  + 

0.56 

Rockville  Dam* 

3  April  1966 

Mine  Shaft 

4.51  + 

0.93 

Tan 

3  June  1966 

Tuff 

5.56  + 

0.49 

Yuba 

5  June  1963 

Tuff 

4.36  + 

0.49 

•Chemical  Explosion 

Table  3,  Data  From  Explosions  Analyzed 
by  the  Seismic  Data  Laboratory. 


A.  General 

The  basic  aim  of  the  VELA  UNIFORM  Program  has  been, 
and  is,  "to  obtain  at  the  earliest  practical  date  a  system 
for  the  detection  of  nuclear  explosions  underground..."  im¬ 
plicit  in  this  objective  is  the  identification  of  a  given 
seismic  event  at  a  nuclear  explosion  once  it  has  been  detected. 
Thus,  the  basic  problem  restated  is  twofold  -  detection  and 
identification.  Tc  the  basic  problem  we  must  alto  a^u  the 
necessity  of  developing  a  system  to  implement  ox,  i  solution. 

Research  progress  in  the  Seismic  Data  Laboratory  has 
paralleled  that  in  the  seismological  community  and  has  ’  ed  in 
most  of  the  specific  areas  subject  to  SDL  attention.  Early 
research  was  restricted  to  recordings  of  regional  and  near 
regional  events.  As  progress  was  made  in  the  solution  of  the 
detection  and  identification  problem  at  these  distances,  atten¬ 
tion  shifted  to  recordings  of  explosion  and  earthquake  seismo¬ 
grams  at  teleseismic  distances. 

During  this  reporting  period,  emphasis  has  been  placed 
on  detection  and  identification,  utilizing  data  recorded  at 
more  sophisticated  facilities,  such  as  the  LASA  in  Montana, 

and  utilizing  more  sophisticated  techniques  for  analyzing  such 
data. 

Following  is  a  discussion  of  the  research  undertaken 
at  the  SDL  for  which  results  were  reported. 

B .  LASA  Studies 

1.  LASA  Signal  and  Noise  Amplitudes 

The  objectives  of  this  study  were  tc  determine  the 
variability  of  signal  and  noise  amplitudes  within  individual  sub¬ 
arrays,  and  across  the  total  LASA,  in  addition  to  signal-to-noise 
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ratio  gains  produced  both  by  phased  and  unphasad  subarray  sums. 
Computations  were  made  of  signal  amplitudes,  rms  noise  levels, 
signal-to-noise  ratios  and  phased  and  unphased  s  <barray  summations 
from  Montana  LASA  recordings  of  three  teleseismic  events. 

The  1,325  digital  seismograms  used  in  this  study 
were  recorded  on  six  multiplexed  tapes  (556  Lpi) ;  three  of  these 
tapes  contained  one  cps  sine  ave  calibration  data,  and  the  re¬ 
maining  three  tapes  contained  event  recordings.  Portions  of  the 
event  seismograms,  amounting  to  120  seconds  of  data  (70  seconds 
before  P  arrival  plus  50  seconds  of  signature) ,  were  placed  in 
SDL  library  format;  these  library  tapes,  two  for  each  event,  were 
then  used  as  input  to  the  LASA  amplitude  program. 

Following  the  data  formatting  stage,  the  seismograms 
were  bandpass  filtered  and  ^magnified  (corrected  for  system  magni¬ 
fication  at  one  cps)  . 

Figure  6  presents  the  pre-filtered  outputs  of  the 
subarray  center  seismometer  which  are  typical  (in  waveform)  of 
those  recorded  throughout  LASA. 

Figure:  7  and  8  Show  subarray  amplitude  and  S/N 
variations,  respectively,  computed  from  the  04  November  1965  data. 
The  mean  signal  amplitude  of  13  subarrays  is  eight  millimicrons, 
whereas  subarray  E3  has  a  mean  amplitude  of  about  20  mjj  (about 
3  times  the  LASA  mean).  Corresponding  S/N  (Figure  8)  also  show 
that  subarruy  E3  is  greater  than  the  mean  by  a  factor  of  three. 

The  subarray  standard  errors  for  signal  amplitudes 
for  the  three  events  used  averaged  about  0.2,  signifying  that 
events  producing  large  mean  signal  amplitudes  have  correspondingly 
large  standard  deviations. 
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Figure  6.  Typical  LASA  Seismograms  for  4  November  1965  - 

Argentina  Event 
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Figure  7.  Subarray  Amplitude  Variation*  for  4  November  1965  Event 
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Figure  8.  Subarray  S/N  Variations  for  4  November  1965  Event 


From  the  foregoing  diecueeion  of  LAS A  eignal  end 
noise  amplitude*  for  three  earthquake*,  the  following  i*  concluded: 

a.  Signal  amplitude*  vary  within  a  subarray  by 
a  factor  of  three  or  four; 

b.  Subarray  mean  aignal  amplitude*  vary  across 
LASA  by  a  factor  of  three  or  four; 

c.  Noise  levels  (rm»)  across  the  entire  1ASA  appear 
to  be  fairly  uniform  for  a  given  events 

d.  Although  the  mean  signal  and  noise  arr^litudes 
computed  from  seismograms  recorded  at  the  sub¬ 
array  centers  are  lower  than  their  corresponding 
values  computed  from  all  LASA  seismograms,  the 
S/N  values  are  about  equal; 

o.  There  are  no  apparent  similarities  among  the 

observed  amplitude  variations,  e.g.,  the  maximum 
amplitudes  computed  from  the  three  events  are 
recorded  by  different  subarrays,  although  two 
of  the  events  are  from  the  same  direction  and 
within  1000  kilometer*  of  on*  another; 

f.  Phased  summations  of  the  subarray  signal  am¬ 
plitudes  are  about  2  db  better  than  corresponding 
unphased  sums,  but  they  are  slightly  smaller  than 
the  subarray  mean  amplitudes; 

g.  Phased  summations  of  the  subarray  noise  levels 
(rros)  are  about  the  same  (actually  slightly 
smaller)  as  the  corresponding  unphased  sums, 
but  they  are  about  6  db  lower  than  the  subarray 
mean  noise  levels;  and 

h.  Subarrays  S/N  from  the  phased  summations  are 
about  3  db  better  than  the  unphased  sums,  and 
they  are  about  6  db  better  than  the  subarray 
mean  S/N, 


2.  Maximum  Likelihood  Filtering  of  LASA  Seismograms 

The  objective  of  this  study  was  to  compare  the 
effectiveness  and  practicality  of  the  maximum-likelihood  filtering 
technique  with  straightforward  time-shifted  array  summation,  i.e., 
beamforming.  We  were  interested  specifically  in  signal-to-noise 
ratio  improvement,  distortion  of  signal  waveform,  and  cost  in 
digital  computing  time. 

The  data  were  short-period  LASA  seismograms  contain¬ 
ing  P  arrivals  from  two  teleseismic  events,  10  November  1965, 

Al'  itian  Islands,  and  25  November  1965,  Kamchatka.  Realizable 
and  symmetrical  maximum-likelihood  filters,  two  seconds  long  and 
based  on  fitting  ii.tervals  consisting  of  150  seconds  of  bandpass 
filtered  (0.4  -  3.0  cps)  data,  were  used  in  processing  the  data. 

The  results  of  the  analysis,  summarized  in  Table  4 
and  illustrated  by  Figures  9  and  10,  indicate  that  within  a  sub- 
array  the  maximum-likelihood  filter  gives  about  1  db  moro  S/N 
improvement  than  the  phased  sum,  when  the  noise  reference  is 
outside  the  fitting  interval.  The  symmetric  maximum-likelihood 
filter  performs  about  2  db  better  than  its  realizable  counterpart. 

Moreover,  when  individual  subarrays  were  processed  by 
the  realizable  filters,  and  the  21  resulting  channels  were  summed, 
the  S/N  improvement  was  2-3  db  over  simply  phasing  and  summing  the 
individual  subarray  sums;  this  improvement  referred  to  noise  out¬ 
side  the  fitting  interval. 

The  symmetric  maximum* likelihood  filter  degrades 
output  signal  waveform  when  signals  on  the  input  channels  arc  not 
identical.  The  degradation,  which  is  usually  not  serious  within 
a  wubarray,  takes  the  form  either  of  precursors  to  the  signal  or 
diminution  of  the  first  motion. 
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Table  4.  Susiiary  of  Signal/feioiae  Ratio  Iaprovesaent 
(decibela  above  bandpaaa-filtered  a  ingle  aeiaaoseter) 


i 


3.  Maximum- Likelihood  Filtering  of  LASA  Noise  Seismograms 

The  purpoae  of  thia  atudy  wai  to  determine  the  effect- 
ivenea*  of  realizable  (MLR)  and  symmetrical  (HLS)  maximum- likelihood 
filters,  compared  with  direct  aummation,  in  reducing  the  rme  noiae 
level  inaide  of,  immediately  adjacent  to,  and  aeveral  roinutee  away 
from  the  fitting  interval:  moreover,  we  were  intereated  in  the 
comparative  behavior  of  maximum-likelihood  filtere  computed  from 
variable  lengtha  of  both  raw  and  bandpaea  filtered  inpute. 

Realizable  and  symmetrical  maximum-likelihood  filters, 
tvo  seconds  long  and  baaed  on  fitting  intervale  coneieting  of  75. 

150,  180,  and  300  seconds  of  both  raw  and  bandpass  filtered  inputs, 
were  computed  using  seismic  noiae  recorded  during  the  night  of  25 
March  1966,  by  nineteen  sensors  in  the  AO  aubarray  at  the  Montara 

LAS  A. 

Results  showed  that,  excluding  the  realizable  filter 
based  on  raw  data,  the  effect  of  maximum-likelihood  filtering  noise 
samples  lying  outside  the  fitting  interval  is  to  produce  noise 
levels  1  to  3  db  lower  than  the  corresponding  direct  sum  of  band¬ 
pass  filtered  data.  Moreover,  the  filters  lower  the  rms  noise 
level  slightly  more  than  the  square  root  of  N  in  the  data  sample 
lying  outside  of,  but  immediately  adjacent  to,  the  fitting  interval, 
but  fail  to  yield  the  square  root  oiT  N  noise  reduction  if  applied 
to  noise  fields  which  are  soveral  minutes  removed  from  the  interval. 

Depending  upon  the  nature  cf  the  input  data  (raw  or 
pre-f iltered) ,  the  symmetrical  filter  is  1  to  5  db  better  than  the 
corresponding  realizable  filter,  but  its  output  appears  to  be 
approximately  100  degrees  out-of-phase  with  its  realizable  counter¬ 
part. 
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It  was  alao  found  that  the  fitting  interval  is  not  a 
critical  parameter  in  the  deaign  of  maximum-likelihood  filters  for 
lengths  equal  to  or  greater  than  150  seconds. 

A.  The  Crosswise  Sum  Method  for  Maximum-Likelihood 

Filtering  of  LASA  Seismograms 

In  previous  maximum-likelihood  studies  of  data  from 
the  Montana  LASA,  individual  subarrays  were  processed  and  then 
these  subarray  outputs  were  combined  to  process  all  of  LASA.  It 
may  be  advantageous  to  phase  and  sum,  across  LASA,  one  output 
from  each  subairay  (say,  all  those  sensors  having  the  same  iden¬ 
tification  number  —  e.g.,  the  phased  sum  of  all  the  10‘s,  tho 
phased  sum  of  all  the  21 ’a,  etc.)  and  then  to  use  the  resulting  25 
traces  as  input  to  a  maximum- likelihood  filtering  process. 

Adding  together  sample i  of  noise  from  all  subarrays 
by  combining  outputs  from  corresponoing  elements  from  each  subarray, 
would  produce  25  traces  with  more  noise  correlation  than  the  21 
subarray  outputs  themselves.  Moreover,  adding  together  samples 
from  each  subarray  would  produce  25  signals  of  similar  waveform. 
Finally,  because  only  one  maximum  likelihood  filter  would  have  to 
be  designed,  the  crosswise  sum  method  would  save  a  considerable 
amount  of  machine  time. 

Accordingly,  it  was  decided  to  process  the  two 
event*  reported  earlier  (10  November  1965  and  25  November  1965) 
by  this  crosswise  sum  method,  and  to  make  the  necessary  compari¬ 
sons  between  results. 
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Results  showed  that  the  crosswise  sum  method  produces 
results  within  approximately  2  db  of  the  maximum  gain  indicated 
by  the  previous  analysis.  However,  whereas  in  previous  results 
the  subarray  maximum-likelihood  filters  produced  gajfns  of  0  to  1.5  db 
over  the  corresponding  phased  sums  for  noise  measured  outside  the 
fitting  interval,  the  crosswise  sum  maximum-likelihood  filter  seems 
to  lose  from  1  to  3  db  relative  to  this  corresponding  phased  sum. 

This  is  in  accord  with  previous  results  for  maximum-likelihood  fil¬ 

ter*  used  to  combine  subarray  outputs,  which  lost  from  1.5  to  4  db 
relative  to  their  corresponding  phased  sums  of  noise  measured  out¬ 
side  the  fitting  interval. 

The  signals  on  the  traces  input  to  the  crosswise  sum 
method  are  remarkably  similar  in  waveform,  and  they  are  egual  in 
amplitude.  This  is  reflected  in  the  preservation  of  signal  ampli¬ 
tude  through  the  maximum-likelihood  filter.  The  symmetrical  filter 
does  not  introduce  a  precursor.  In  both  realizable  and  symmetrical 
cases,  the  signal  appears  as  it  does  on  the  corresponding  phased  sum. 

For  each  of  the  inputs  to  the  crosswise  sum  raaximum- 

likelihced  filters,  a  theoretical  gain  of  13.2  db  is  predicted  if 
the  noise  is  completely  uncorrelated,  because  here  N  «  21  and  the 
square  root  of  N  improvement  is  13.2  db.  For  the  25  November  event, 
this  is  achieved  for  the  average  of  the  input  traces.  For  the 
10  November  event,  it  is  lower  by  3  db. 

5.  Relative  Travel-Time  Anomalies  at  LASA  and  the 

Location  of  Epicenters  Using  "SHIFT" 

The  entire  set  of  1B0  events  selected  for  tie  measure¬ 
ment  of  anomalies  were  recorded  with  good  signal-to-noise  ratios  on 
the  films.  If  any  of  the  subarrays  recorded  signals  were  difficult 
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to  read,  they  were  eliminated  from  all  computations.  However, 
if  tne  signals  were  strong,  and  if  confidence  was  placed  in  the 
readings,  the  calculated  anomalies  were  keot  regardless  of  their 

agreement  with  other  results. 

Whenever  possible,  the  readings  were  the  arrival 

times  of  first  motion  of  P.  However,  the  number  of  events  re¬ 
corded  with  sufficiently  distinct  onsets  was  small  and  the  data 
accumulated  too  slowly.  Because  all  of  the  instruments  at  LASA 
are  the  sam» .  first  peaks  c.  troughs  of  a  number  of  events  were 
also  read  to  increase  the  rat  of  accumulation.  When  the  anomalies 
obtained  by  using  first  motions  were  closely  compared  with  those 
obtained  using  first  extrema,  no  consistent  differences  were  ap¬ 
parent.  This  result  implies  that  any  differences  in  the  phase 
response  characteristics  between  the  subarray  center  instruments 
are  within  the  precision  of  the  measurements. 

The  event  locations,  necessary  for  computing  ex¬ 
pected  travel-times,  were  those  listed  on  the  "Preliminary  Deter¬ 
mination  of  Epicenters"  ( PDE)  cards  supplied  by  the  U.S.  Depart¬ 
ment  of  Commerce,  Environmental  Science  Services  AdministraJ ion 
(U.S.  Co?st  and  Geodetic  Survey).  These  locations  are  routinely 
reported  to  0.1°  latitude  and  longitude.  Merely  for  convenience 
the  center  subarray  of  LASA,  AO,  was  selected  as  the  reference 

station  for  all  calculations. 

In  general,  travel-time  anomalies  are  dependent 

upon  both  event  direction  and  distance.  It  is  necessary  to 
measure  anomalies  from  a  large  number  of  events  occurri.x  in  all 
regions  of  possible  interest.  More  consistency  is  obtained  if 
the  results  are  then  separated  into  event  direction  of  approach 
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and  plotted  as  a  function  of  distance.  It  was  seen  that  if  the 
anomalies  were  plotted  as  a  function  of  direction  with  no  regard 
for  distance,  the  results  would  be  quite  scattered.  Similarly, 
if  they  were  plotted  only  as  a  function  of  distance,  the  same 
large  scatter  would  result.  Both  distance  and  direction  are 

important  and  must  be  separated. 

Although  the  anomalies  are  plotted  by  direction  and 

distance,  they  are  more  precisely  plotted  by  region.  The  term 
"region"  as  used  here  is  defined  as  that  limited  area  from  which 
all  events  will  yield  the  same  anomalies. 

The  travel-time  anomalies  at  LASA  were  seen  to  be 
large,  even  for  stations  quite  near  the  reference,  AO.  For  example, 
subarray  D4  recordec  signals  0.7  seconds  late  with  respect  to  AO 
for  events  southeast  of  LASA  at  a  distance  of  about  8000  km.  The 
physical  causes  of  the  anomalies  are  presently  unknown,  but  it  is 
probable  that  differing  geological  conditions  in  the  crust  and  the 
upper  mantle  beneath  each  of  the  subarrays  is  responsible. 

The  results  showed  scatter  within  some  regions  much 
larger  than  would  be  expected  from  reading  errors  alone.  The  only 
reasonable  cause  of  the  observed  anomaly  variation  within  a  region 
is  an  error  in  the  assumed  location  of  the  event.  In  this  study, 
all  of  the  epicenters  and  depths  used  were  those  listed  on  the  PDt 
cards.  These  locations  are  meant  to  be  preliminary,  and  they  are 
accepted  to  be  as  such.  The  reported  latitudes  and  longitudes 
"have  an  estimated  accuracy  of  a  few  tenths  of  a  degree  and  25  km 
in  depth"  (U.S.C.&.G.S  first  PDE  card).  An  error  of  a  "few  tenths 
of  a  degree" can  easily  explain  most  of  the  anomaly  variations. 
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The  next  logical  question  ies  Can  the  original  epi¬ 
centers  within  a  region  be  shifted  by  some  acceptable  distance  such 
that  the  observed  scatter  about  a  mean  is  reduced?  To  aid  in  an¬ 
swering  this,  a  computer  program, SHIFT, was  written  which  reduces 
the  anomaly  errors  by  shifting  the  locations.  An  example  of  the 
anomaly  improvement  at  just  one  of  the  20  subarrays  using  Kurile 
Islands  events  is  shown  in  Figure  11,  where  the  results  in  the 
upper  part  indicate  the  computed  anomalies  from  the  original  epi¬ 
centers  and  the  results  in  the  lower  part  those  after  the  shifts. 

6.  LAS A  Travel-Time  Data  at  the  SDL 

In  connection  with  its  many  LAS A  analyses,  the  SDL 
is  accumulating  relative  travel-times  at  the  Montana  Large  Aperture 
Seismic  Array  (LASA)  with  a  view  towards  computing  travel-time 
anomalies  at  the  21  subarrays.  In  the  belief  that  these  raw  data 
are  in  demand  and  may  be  of  use  to  the  seismic  community  in  general 
a  report  was  issued  which  contains  LASA  travel-time  data  for  approx 
imatelv  400  events,  as  read  from  LASA  films. 

Figure  12  shows  the  type  of  data  which  has  been  com¬ 
puted  and  the  following  numbered  explanation,  corresponding  to  the 
circled  number  on  the  figure,  describes  the  method  of  presentation: 

1.  Arbitrary  region  name 

2.  Direction  of  approach 

3.  Distance  range  for  events  included  in  region 

4.  Azimuth  range  for  events  included  in  region 

5.  Evunt  date  and  name 

6.  PDE  latitude 

7.  PDE  longitude 

8.  PDE  depth 

9.  PDE  origin  time 

10.  Arrival  time  in  seconds  (an  arrival  time  of  "0" 
indicated  no  reading  made  at  that  subarray) 

11.  Arrival  time,  hour,  and  minute  at  LASA 
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Figure  12.  Type  of  LASA  Travel-Time  Data  Computed  by  the  SDL 
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7.  LASA  Travel-Time  Anomaliea  for  Various  Epicentral 
Regions 

Approximately  350  telepeisms  were  used  to  compute 
P-wave  travel-time  anomalies,  relative  to  subarray  AO,  at  all 
subarray  center  instruments  at  the  Large  Aperture  Seismic  Array 
(LASA)  in  Montana.  As  the  anomalies  were  dependent  upon  event 
distance  and  azimuth  (hence,  epicentral  region) ,  the  results  were 
separated  into  groups  such  that  the  subarrays  had  anomalies  which 
were  consistent  within  each  defined  region.  The  total  number  of 
regions  was  39. 


The  results  were  arranged  first  by  general  direction 
(beginning  with  the  northwest  and  going  clockwise)  and  second  by 
increasing  epicentral  distance  within  each  directional  group,  as 


follows: 

No. 

Azimuth 

No. 

of 

Distance 

Range 

of 

Direction 

Regions 

Ranqe,  km 

(Deqrees) 

Events 

Northwest 

14 

2000-10800 

292-320 

155 

North 

2 

9300-10900 

349-003 

14 

North-nor  thea  st 

1 

9300-  9900 

018-025 

3 

Northeast 

2 

4800-  9800 

032-054 

12 

East 

1 

6000-  6100 

070-071 

2 

East-southeast 

2 

4800-  5800 

112-120 

6 

Southeast 

9 

3400-  9800 

133-165 

90 

South 

3 

2800-  9100 

181-188 

12 

West- southwest 

2 

9500-10800 

239-246 

34 

Undefined 

3 

Atlantic  Ridge 

5000-  9900 

083-101 

6 

Continental  U. S. 

800-  1700 

124-280 

5 

Miscellaneous 

5300-10000 

328-178 

6 
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The  following  observations  are  made  concerning 
travel-time  anomalies  at  LASA: 

1.  The  anomaly  variations  between  regions  measured 
at  a  single  subarray  are  as  high  as  1.43  seconds.  For  example, 
subarray  F4  has  an  average  anomaly  of  +0.65  sec  (N  =  8,  n  ^  ■  0.06) 
for  the  Eastern  Mexico  region  and  an  average  of  -0.78  sec  (N  **  17, 

■=  0.05)  for  the  Fiji  Islands  region. 

2.  Subarrays  which  are  quite  near  the  reference 

subarray  A0  can  have  unusually  large  anomalies.  For  example, 
subarray  D4  is  located  30.75  km  from  A0  (center  instrument  to 
center  instrument)  and  has  an  anomaly  relative  to  A0  of  +0.71  sec 
(N  ■  12,  ■  0.06)  for  the  Peru-Bolivia  region. 

3.  The  center  instrument  at  subarray  B2  is  only 

7.50  km  from  the  center  at  A0,  but  it  has  an  anomaly  of  -0.26  sec 

(N  »  11,  a  *  0.07)  for  the  Fiji  Islands  region.  This  result 
B2 

suggests  that  the  time  anomalies  within  one  subarray  (7  km  diameter) 
may  vary  from  element  to  element  perhaps  as  much  as  0.30  sec.  Sig¬ 
nals  within  a  subarray  could  be  significantly  misaligned  if  such 
large  anomalies  existed. 

4.  The  anomalies  are  not  slowly  varying  functions 

of  either  distanco  or  azimuth.  For  example,  subarray  El  has  an 

anomaly  of  +0.78  sec  (N  *  12,  ■  0.06)  computed  from  events 

approaching  from  a  southeasterly  direction  at  7700  km  distance 

(Peru-Bolivia  region),  but  at  9300  km  (Central  Chile-Argentina 

Border  region) ,  whereas  for  events  bearing  about  182°  at  the  same 

distance  (South  and  West  of  Easter  Island)  the  anomaly  is  -0.57 

sec  (N  »  6, o  »  0.09).  Hence  the  anomaly  at  El  changes  by  0.53 

F3 

sec  in  a  distance  of  1600  km,  and  at  F3  it  changes  by  0.68  sec  in 
38°  of  azimuth. 
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5.  The  maximum  anomaly  range  observed  at  LASA  is 
1.72  sec;  average  anomaly  at  subarray  El  is  +0.87  sec  for  events 
occurring  in  the  North  Colombia  region;  average  anomaly  at  subarray 
F3  is  -0.85  sec  for  events  from  the  North  Atlantic  region. 

6.  The  maximum  anomaly  range  for  one  particular 
region  (North  Atlantic)  is  1.18  sec,  where  the  F3  anomaly  is  -0.85 
sec  and  the  Fl  anomaly  is  +0.33  sec. 

7.  In  general,  the  aubarray  anomalies  computed  for 
a  region  are  consistent  within  the  precision  of  the  measurements. 
There  are,  however,  discrepancies.  For  example,  events  occurring 
in  the  North  Colombia  region  yield  anomalies  which  differ  more 
than  expected  from  one  event  to  another,  and  the  subarray  and 
event  standard  deviations  (  c .  and  o  .  respectively)  are  much 
larger  than  those  computed  for  other  regions.  The  probable  reason 
for  these  discrepancies  is  epicenter  mislocation;  i.e.,  these  re¬ 
gions  may  not  have  adequate  seismic  control  such  that  spurious 
locations  result. 

8.  It  is  estimated  that  the  36  regions  (plus  3  un¬ 
defined)  given  in  this  report  yield  adequate  anomaly  information 
for  about  20%  of  all  possible  regions  and  for  about  90%  of  all 
earthquake  regions  within  105°  of  LASA, 

C.  Teleseismic  Signal  Measurements  at  Vertical  Arrays 

A  vertical  array  records  simultaneously  on  several 
transducers  stacked  in  a  deepwell  3  km  or  more  in  depth.  The 
purpose  of  this  study  was  to  investigate  the  possibility  of  re¬ 
ducing  near  surface  reverberations  by  stacking  procedures  similar 
to  those  used  in  reflection  seismology.  By  reciprocity,  the  same 
algorithms  developed  for  stacked  sources  can  be  applied  to  stacked 
receivers  by  reversing  the  sign  of  the  reflection  coefficients. 


a a  follows! 


The  two  deghosting  techniques  applied  are  described 


1.  When  a  surface  or  near  surface  trace  is  avail¬ 
able,  the  deghosted  trace  can  be  constructed  by  simply  shifting  the 
surface  trace  by  K,  multiplying  a  »s  a  ,  and  subtracting  this  re¬ 
sultant  from  either  the  d&epwell  trace  (phased  vertical  array)  or 
the  deepwell  trace  shifted  by  K/2  (non-phased  vertical  array). 

2.  It  would  seem  more  plausible  to  develop  a  de- 
ghosting  process  that  removes  the  echo  without  the  use  of  a  surface 
seismogram,  a  practical  technique  that  would  figuratively  push  the 

ghost  off  the  end  of  the  deepwell  trace. 

Given  the  reflection  coefficient  ft  ,  and  the  deep- 
well  trace  y±  1  the  phased  deghosted  trace  is  formed  as  follows: 
First  a  new  crace  Z ^  ^  is  defined  to  be 

Zi,l  ”  Yi,l  “XU  -  NTx)  '  «  1, 

where  NT  is  equal  to  the  product  of  the  sampling  rate  and  the 
echo  time  delay  (K)  .  This  above  equation  simply  shifts  and  in¬ 
verts  the  ghost  to  a  point  farther  down  the  seismic  record  than 

it  originally  was.  This  iterative  process  is  continued  until  the 
ghost  reflection  is  pushed  off  the  seismic  record. 

The  individual  deghosted  traces  from  the  vertical 
array  are  then  stacked  by  summing.  The  signals  should  now  re¬ 
inforce  since  after  deghosting  they  tend  to  agree  in  waveform 
and  arrival  tine.  The  noise  background  and  oscillations  created 
on  individual  channels  by  the  deghosting  process  will  decrease 
relative  to  the  signal  upon  summation. 

-32- 


4 


In  addition  to  the  linear  summation,  a  correlation 
trace  ia  formed  which  emphasizes  the  in-phaae  component  or  all  n 
channel*  in  the  vertical  array  tracea.  The  correlation  trace  ia 
formed  by  firat  computing  a  running  zero  lag  correlation,  y,,(t) , 
between  the  first  two  tracea,  x^(t)  and  *****  vort*ca^  ®*ray. 

t+T 

y2(t)  ■  2t"  J  xl(u)  X2(u)  dU 

t-T 

Then  the  aubaequent  tracea  in  the  vertical  array  are  correlated  in 

a  airailar  way  with  the  correlation  trace. 

t+T 

y3*t)  "  2T~  I  y3<U)  X3(U)  dU 
t-T 

*  t+T 

yk(t)-irJ  Vi(u>  Vu)  du 

t-T 

For  k  tracea  in  the  vertical  array,  yk(t)  repreaenta  the  final 
correlation  trace  which  ia  displayed. 

Figure  13  shows  the  original  recordings  for  an 
Aleutian  earthquake  recorded  at  the  AP-OK  vertical  array.  The 
aum  trace  and  the  correlation  trace  are  computed  with  no  de¬ 
ghosting  applied  to  the  individual  traces. 

Figure  14  shows  this  same  event  with  the  individual 
traces  deghosted  by  method  I.  Figure  15  shows  this  same  event 
with  the  individual  traces  deghosted  by  method  II. 
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Vertical  Component  Vertical  Array  Measurements  --  Aleutians# 
USCGS  Focus  67  km,  Magnitude  5.2,  Date  6/17/65,  Time  19*05*9,1 


Figure  13 
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Figure  15.  Ceghosted  Vertical  Array  Measurements,  Method  II 

(Reflection  Coefficient  0.7) 


Results  shew  that  vertical  array  processing,  based 
on  previous  experience  with  stacked  sources,  reduces  the  rever¬ 
berations  and  complexity  of  the  coda  caused  by  near  surface  re¬ 
verberations  at  the  receiver.  Both  deghosting  methods  performed 
adequately.  The  correlation  trace  for  the  inverse  operation  had 
a  slightly  cleaner  coda  and  this  method  has  t'.«  advantage  of  not 
requiring  a  quiet  surface  trace  and  exact  control  of  instrument 
gain.  The  correlation  traces  used  for  detecting  the  upgoing  P- 
pulses  impose  the  strongest  possible  requirement  that  the  signal 
be  fixed  jointly  on  all  channels. 

D.  Rayleigh  Wave  Rejection  by  Optimum  Filtering  of 

Vertical  Arrays 

Vertical  arrays  offer  some  intriguing  advantages  over 
other  arrangement  of  seismometers.  Assuming  the  noise  to  be 
composed  of  Rayleigh  waves  and  perhaps  some  mantle-propagating 
P-waves,  and  assuming  the  signals  to  be  teleseismic  P-waves,  the 
noise  and  signals  will  be  recorded  by  a  vertical  array  with  some 
rather  special  characteristics. 

We  can  imagine  that  each  mode  of  Rayleigh  wave  noise 
possesses  a  random  character  common  to  noise  functions  in  general. 
However,  tor  each  Rayleigh  wave  mode  the  response  versus  depth 
and  frequency  relative  to  its  response  on  the  surface  is  predict¬ 
able.  Because  the  relative  depth  variation  is  predictable,  the 
vertical  array  can  be  summed  to  cancel  the  Rayleigh  modes.  In 
contrast,  the  noise  over  a  surface  array  is  more  unpredictable 
and  more  variable  and,  therefore,  more  difficult  to  cancel  by 
array  summation. 


The  purpose  of  this  analysis  was  <o  test  the  theory  of 
applying  maximum  likelihood  filters  in  vertic  1  arrays  in  order 
to  provide  undistorted  estimates  of  the  signal  or  the  various 
Rayleigh  modes  with  the  other  modes  cancelled  out. 

A  vertical  array  of  five  vertical  component  seismometers 
which  existed  in  a  well  at  the  Uinta  Basin  Seismological  Observa¬ 
tory  (UBSO)  was  used  for  testing  this  theory.  The  Rayleigh  dis¬ 
persion  curves  for  this  well  indicated  that  the  fundamental,  first, 
and  second  higher  Rayleigh  modes  would  all  be  supported  over  the 
signal  frequency  rangf  of  0.5  to  2.0  cps.  Since  we  wanted  more 
recording  levels  than  modes,  these  three  Rayleigh  modes  plus  a 
signal  mode  were  the  only  ones  considered  to  exist  in  our  syn¬ 
thetic  modeling  of  this  well. 

Che  optimum  filter,  G,  solutions  are  the  least  squares 
inverse  to  the  Rayleigh  signal  filter  matrix,  H.  Thus,  the  matrix 
products  of  G  and  H  gives  the  identity  matrix.  In  our  test  model 
we  assumed  the  signal  existed  with  the  same  size  and  amplitude  on 
all  traces.  We  derived  optimum  filter  solutions  only  over  the 
signal  range  of  0.5  to  2.0  cps.  For  the  frequencies  below  0.5, 
we  smoothed  the  frequency  responses  to  zero  according  to  a  sinusoidal 
gain  function  for  both  the  H  and  G  frequency  responses.  We  smoothed 
these  responses  to  zero  in  a  similar  way  over  the  range  from  2  to 
5  cps.  Also,  the  orthogonality  of  our  optimum  filter  solutions 
were  maintained  over  the  entire  frequency  range. 

In  order  to  demonstrate  that  our  optimum  filter  solutions 
were  correct,  we  created  some  synthetic  data  which  conformed  to 
the  well  log  analysis  we  assumed  for  the  UBO  well.  These  Rayleigh 
mode  traces  were  combined  with  an  artificial  signal  to  produce  a 
mixture  of  signal  and  Rayleigh  modes  expected  from  this  well. 
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Figure  lb  shows  the  data  from  the  five  levels  of  the 
vertical  array  if  the  signal  plus  the  fundamental  and  the  first 
and  second  higher  Rayleigh  modes  have  been  added  into  the  vertical 
array  data.  The  outputs  for  the  signal  and  the  three  Rayleigh 
modes  a*-e  reproduced  in  size  and  waveform  within  a  few  percent  of 
the  expected  results.  The  signal  trace  shows  better  than  a  20  db 
improvement  over  the  best  signal/noise  ratio  available  in  the 
vertical  array.  Thus,  the  optimum  filters  perform  as  expected 
when  the  data  properties  match  the  dispersion  analysis  of  the  well 
^  log. 

These  solutions  were  then  applied  to  actual  noise  data 
recorded  at  the  vertical  array.  All  of  the  optimum  filters  used 
reguired  inputs  from  all  five  levels  of  noise  data.  However,  we 
could  ask  for  estimates  of  only  a  single  mode  or  for  estimates 
of  several  modes  in  the  filter  outputs.  The  more  modes  a  set  of 
filters  estimate,  the  larger  the  filter  gains  will  be.  It  follows 
that  any  errors  resulting  from  incorrect  dispersion  analysis,  in¬ 
strument  calibrations,  or  non-parallel  layering  will  be  amplified 
much  more  in  these  filters  with  fewer  degrees  of  freedom.  Results 
showed  that  the  signal  estimate  outputs  gave  progressively  in¬ 
creasing  low  frequency  errors  on  the  estimates  of  all  modes  as 
the  degrees  of  freedom  decreased  (i.e.,  as  the  number  of  outputs 
increased) . 

It  was  concluded  that  our  optimum  filter  solutions  have 
the  following  features: 

1.  These  filters  are  indeed  orthogonal  and  separate 
the  noise  and  signal  modes  as  planned. 

2.  These  optimum  filters  extend  for  400  points  in  time 
since  the  frequency  interval  chosen  was  0.5  cps. 
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Optimum  Filters  Estimating  Four  Nodes 
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3.  The  gain  on  the  synthetic  examples  was  greater 
than  20  db  with  solutions  restricted  to  the  0.5  -  2.0  frequency 
range.  More  resolution  in  frequency  could  increase  this  figure. 

4.  The  optimum  filters  are  zero  phase  shift  filters. 
Therefore,  they  are  algebraically  additive. 

5.  The  optimum  filters  become  larger  in  gains  (both 
positive  and  negative)  as  degrees  of  freedom  decrease. 

6.  The  optimum  filters  become  larger  in  gains  (both 
positive  and  negative)  as  the  aperture  of  array  goes  down.  Thuc 
for  low  frequencies  the  solutions  tend  to  become  unstable. 

7.  Extra  modes  not  considered  and  errors  in  assumptions 
cause  errors  in  the  output.  Errors  in  assumptions  can  include  an 
incorrect  well  log,  non-parallel  layering  in  surrounding  medium, 
and  incorrect  calibrations  of  seismometers. 

8.  Extra  degrees  of  freedom  are  needed  to  cut  gains  of 

optimum  filters  and  make  optimum  solutions  more  tolerant  of  errors. 

9.  Extra  degrees  of  freedom  are  best  obtained  by  in¬ 
creasing  the  number  of  seismometers  in  the  vertical  array. 

10.  More  stable  solutions  (i.e.,  optimum  filters  with 
lower  gain)  will  be  obtained  from  the  deeper  vertical  arrays  which 
have  the  seismometers  distributed  rather  uniformly  throughout  the 
array. 

E.  Detection  of  Surface  Waves  at  Teleseismic  Distances 

We  present  here  a  matched  filter  approach  for  distin¬ 
guishing  weak  teleseismic  surface  wave  signals  from  background 
noise.  The  method  discriminates  against  events  not  located  in  a 
particular  source  region  of  interest  and  provides  estimates  of 
magnitude  and  radiation  pattern,  when  a  number  of  recording  sta¬ 
tions  are  available. 
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Basically  the  matched  filter  approach  amounts  simply 
to  searching  a  record  x(t)  for  a  known  waveform  y(t).  In  par¬ 
ticular,  it  is  assumed  that  x(t)  -  ay(t)  +  n(t)  where  a  is  a 
constant  and  n(t)  is  a  random  noise  process.  If  no  further 
assumptions  are  made  with  regard  to  the  nature  of  the  noise 
process,  one  can  determine  the  least-squares  estimate  of  a  which 
minimizes  J(a,T  )  ■  fx(t  +t  )  -  a  y(t)  2  where  the  summation 
is  over  the  length  of  y(t)  .  In  the  test  series  x,  which  ii»  of 
longer  duration  than  y,  the  lag  *  indicates  at  what  point  in  x 

the  comparison  is  begun.  The  value  of  a  obtained  by  setting  £_J  -  0 
is  da 

a  «  E  x(t  +  T)  y ( t) /  T  y2(t)  (1) 

Thus,  the  matched  filter  in  this  case  is  simply  che  waveform  y(t) 
and  the  matched  filter  output  at  lag  r  is 

£x(t  +  j  )  y ( t )  ■  Cxy  (  T  )  (2) 

The  coherency  at  lag  T  is  given  by 

C(  T)  -  £  x(t  +  T)  y ( t) /  [E  x2(t  +  t)  ‘  E  y2(t)J  *  (3) 

It  always  is  bounded  by  -1  <  C(t)  <  1.  The  maximum  in  the  en¬ 
velope  of  C(  t)  occurs  at  the  value  of  T  where  x  and  y  match 
beat  in  the  least-squares  sense.  The  C  (T )  is  the  correlation 

uill  X 

coefficient. 

In  the  present  application  we  have  used  this  simple 
least  squares  approach.  Estimates  of  a,  C  (  t),  and  C  (  t) 
calculated  in  this  way  allow  us  to  study  test  cases  and  unknown 
events  with  a  minimum  of  assumptions. 


The  matched  filter  approach  requies  that  one  knows 
the  desired  or  expected  waveform  in  order  to  search  effectively 
for  that  waveform  in  a  noisy  record.  in  order  to  obtain  a  suit 


able  "expected"  waveform  and  ac  the  same  time  assure  that  propa¬ 
gation  effects  would  be  properly  accounted  for.  we  chose  the 
surface  wave  from  a  larger  event  in  the  source  region  of  interest 
as  the  filter  y(t).  Thus,  no  matter  how  complicated  the  paths 

t  propdgatlon  are'  fche  paths  for  other  events  in  the  same  source 
region  will  be  nearly  coincident  with  those  of  the  larger  -vent 
which  produced  y(t),  so  the  only  major  differences  in  recorded 
waveform  will  be  source  differences.  This  coincidence  of  travel 
path  simply  assures  that  the  earth's  transfer  function  is  the  same 
for  all  events  occurring  in  that  source  region;  it  does  not  mean 
that  dispersion  during  propagation  to  teleseismic  distances  is 
unimportant.  To  the  contrary,  it  is  this  physical  effect  (dis¬ 
persion)  which  is  vital  to  the  success  of  the  technique,  since 
it  transforms  the  source  pulse  into  an  oscillatory  signal  which 
is  of  long  duration.  As  the  duration  of  the  signal  increases, 
random  correlations  of  the  signal  with  the  noise  become  poorer, 
with  the  result  that  the  "false  alarm"  level  is  reduced.  Prom 
this  standpoint  the  larger  the  epicentral  distance  the  better  the 
method  should  work,  since  the  signal  increases  in  duration,  the 
lar  -er  the  distance.  However,  the  energy  density  of  the  signal 
is  reduced  due  to  dispersior  and  attonuation  during  propagation, 
so  that  the  signal  to  noise  ratio  is  reduced  with  increasing 
distance.  Therefore,  for  a  given  transfer  function  and  noise 
level,  one  would  expect  there  to  be  an  optimum  observing  distance 
for  detecting  events  using  the  matched  filter.  But,  because  of 
the  different  noise  levels  among  available  st^'.ons  and  the 
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variety  of  propagation  paths  to  these  stations,  it  does  not 
appear  to  be  practical  to  determine  these  optimum  observing 
distances  experimentally.  What  we  can  estimate,  however,  is 
the  minimum  S/N  at  which  the  technique  will  detect  observed  or 
expected  signal  waveforms. 

Because  of  dispersion  and  frequency-dependent  attenua¬ 
tion  with  distance,  the  signal  waveform  will  in  general  be  dif¬ 
ferent  at  each  etation  and  thus  a  different  filter  for  each  station 
and  source  region  must  be  used.  However,  this  presents  no  problem 
since  a  single  large  event  in  a  given  source  region  provides  for 
each  station  a  filter  suitable  for  detecting  other  events  located 
in  that  region,  as  well  as  the  temporal  pattern  of  arrival  times 
at  the  stations.  This  pattern  and  the  fact  that  the  filters  are 
different  for  each  station  can  be  used  to  advantage  in  discrimina¬ 
ting  against  events  outside  the  source  region  of  interest,  since 
one  or  more  of  the  decision  criteria  listed  earlier  will  not  be 
satisfied;  in  particular,  the  correlation  peak  values  will  be 
degraded,  and  the  arrival  times  of  these  peaks  will  not  produce 
the  appropriate  pattern  of  arrival  times  at  the  observing  stations. 

When  a  small  event  has  been  detected  by  this  method, 
then  the  question  arises  as  to  what  details  we  can  extract  about 
the  source.  The  estimates  of  a  provide  a  convenient  means  of  com¬ 
paring  the  magnitude  of  the  small  event  to  that  of  the  latger  one, 

*2 

since  a  is  an  estimate  of  the  ratio  of  the  energy  in  the  small 
event  to  that  in  y(t).  Some  care  must  be  taken,  however,  if 
the  reference  event  exhibits  a  strong  radiation  pattern  which  is 
different  from  that  of  the  smaller  event.  In  fact,  the  radiation 
pattern  of  the  small  event  relative  to  the  reference  event  can  be 
obtained  by  plotting  a/^  y^(t)  vs.  source-station  azimuth. 
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Comparisons  of  a  vs.  distance  along  a  given  azimuth 
allow  us  to  estimate  roughly  the  excitation  spectrum  of  the 
small  event  compared  to  the  larger  one,  since  the  attenuation 
with  distance  is  frequency-dependent,  which  implies  that  a 
should  vary  with  distance  in  a  way  which  reflects  the  shape  of 
the  spectrum  of  the  small  event  compared  wi'ih  that  of  the  ref¬ 
erence  event.  Only  if  their  spectra  are  identical  in  shape  will 
a  be  invariant  with  distance.  Since  this  attenuation  factor  is 
of  the  form  exp  [_  — g <uj )  A  j,  where  g(u>  )  is  an  increasing  function 
of  frequency  and  ±  is  epicentral  distance,  %  will  d-screase  with 
distance  if  the  spectrum  of  the  small  event  is  peaked  at  a  higher 
frequency  than  the  reference  event. 

Thus  potentially  the  matched  filter  method  permits  us 

to : 

1*  detect  the  surface  waves  from  a  source  region  of 
interest,  while  rejecting  events  outside  the  source  region; 

2.  make  magnitude  estimates  for  small  events  relative 
to  the  reference  event; 

3.  outline  radiation  patterns  for  small  events; 


4.  estimate  spectral  shape  relative  to  the  reference 


event. 


Figures  17,  18,  and  19  show  the  results  of  adding  the 

signal  shown  at  the  top  (an  actual  oceanic  Rayleigh  wave  train 

recorded  at  4686  km  distance)  to  the  noise  trace  at  n  number  of 

arrival  times,  with  the  indicated  signal/noise,  S/N,  and  using 

the  signal  au  the  matched  filter.  We  define  S/N  as  m^xly(t)  J'RMS 

.  The  number  in  parentheses  is  an  alternate  value  of 

signai/noise  (S/N)  given  by  RMS  ~y(t)j  /RMS  £  n(t)  j.  The  values 

of  cl*  a  are  printed  next  to  C  trace,  and  the  correlation  co- 

xy 

e^ficients  CC  next  to  the  C  trace.  The  arrows  indicate  the 
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Figure  17.  Real  Signal  Test  Case 
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Figure  19.  Real  Signal  Teat  Case 


beginning  time  for  the  signal  in  each  instance.  In  this  example 
the  signal  is  detected  for  S/N  (S/N)  as  low  as  .35  (.15),  where 
visual  detection  is  impossible.  Notice  also  that  the  4's  are 
approximately  in  proportion  to  the  S/N  used  as  expected. 

The  reason  the  results  appear  to  be  as  good  or  better 
for  S/N  *  .35  compared  to  S/N  ■  .5  is  that  the  noise  level  was 
not  uniform  over  the  entire  noise  trace  used  to  establish  the 
RMS  of  the  noise  and  the  .35  signal  was  buried  at  a  position 
where  the  actual  noise  level  over  the  filter  length  was  low 
relative  to  that  where  the  .5  signal  was  buried. 

On  the  basis  of  this  investigation  we  conclude  that: 

1.  In  principle  the  matched  filter  approach  can  be 

used: 

(a)  to  detect  weak  surface  wave  signals  for  signal 
to  noise  ratios  as  low  as  about  0.35; 

(b)  to  make  reliable  relative  magnitude  estimates 
for  S/N  as  low  as  .5; 

(c)  to  determine  radiation  pattern  relative  to  a 
reference  event; 

(d)  to  estimate  the  general  shape  of  a  small  event's 
amplitude  spectrum  relative  to  that  of  the  reference  event. 

2.  The  results  presented  for  a  number  of  actual  events 
have  demonstrated  the  practical  applicability  and  usefulness  of 
the  matched  filter  approach  in  studying  weak  teleseisraic  surface 
waves. 


3.  The  matched  filter  used  in  conjunction  with  the 
surface  wave  arrival  time  patterns  appropriate  for  different 
source  regions  of  interest  allows  one  to  discriminate  against 
all  but  those  events  in  the  particular  source  region  of  interest. 
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F.  Rectilinear  Motion  Detection 


The  studies  of  rectilinear  motion  detection  are  summarized 
in  three  sections.  The  first  section  defines  the  rectilinear 
motion  detectors  (REMODE)  as  non-linear  filters  based  upon  a  cor¬ 
relation  function  between  the  vertical  and  radial  components. 

The  section  described  three  variations  of  kEMODE,  exhibits  the 
application  of  these  filters  to  a  suite  of  real  seismograms  from 
a  Peruvian  earthquake,  and  comments  on  the  different  behavior  of 
the  various  types. 

The  second  section  describes  some  controlled  tests  of  an 
actual  seismic  noise.  This  section  shows  for  one  particular 
REMODE  operator  that: 

a.  with  input  S/N  ratios  of  18  db  or  more  tne  output 
signal  size  is  not  changed  by  REMODE  and  an  average  S/N  ratio 
gain  of  9  db  results; 

V.  with  input  S/N  ratios  of  9  db  to  18  db  the  output 
signal  size  rr.ay  be  changed  by  REMODE  but  a  gain  in  S/N  does  occur; 

c.  with  input  S/N  ratios  less  than  9  db  the  output  signal 
size  is  changed  by  REMODE  and  a  loss  in  S/N  occurs. 

The  third  section  describes  some  controlled  tests  of  an 
actual  signal  buried  in  a  portion  of  polarized  noise  generated 
in  the  coda  of  a  large  earthquake.  The  study  compares  REMODE 
operators  with  various  parameter  settings  in  detecting  a  recti¬ 
linear  signal  in  both  shear  and  rectilinearly  polarized  noise. 

The  results  show  how  much  more  effective  REMODE  is  in  detecting 
rectilinearly  polarized  signals  in  shear  noise  than  in  recti¬ 
linearly  polarized  compressions!  noise. 
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1.  Applications  and  Development  of  Polarization 

(REMODE)  Filters 

The  polarization  of  particle  motion  in  three-component 
seismic  recordings  is  useful  in  the  detection  and  identification 
of  seismic  phases.  Operators  have  been  described  which  use  the 
behavior  of  the  integrated  product  of  two  components  to  detect 
and  identify  phases.  More  general  polarization  filters,  called 
REMODE,  have  been  developed  in  which  the  filter  functions  are  based 
on  the  correlation  functions  of  two  components. 

If  the  original  vertical  and  radial  components  z(t) 
and  r(t)  of  a  seismic  recording  are  used  as  inputs  to  a  REMODE 
process,  the  output  traces  z'(t)  and  r’(t)  are  the  input  traces 
filtered  by  the  time-varying  REMODE  filter  $  : 

t  +T/2 


z  1  ( t  ;  WT) 


-  I 


t  -T/2 
o 


*  t  ,w  (t  -t)  z ( t)  dt 
o  o 


where  i  is  a  filter  function  of  length  T,  derived  from  the 

lagged  cross-products  between  z(t)  and  r(t)  in  a  time  window  of 
length  W  centered  on  tQ.  The  filter  function  varies  as  the  window 
is  moved  along  the  seismogram. 

The  filter  functions  differ  as  follows.  The  filter 
functions  for  REMODE  2  are  time-varying  estimates  of  the  real 
parts  of  cross-correlation  functions. 

The  filter  functions  for  REMODE  3  are  the  REMODE  2 
normalized  by  time-varying  estimates  of  average  power  in 
the  input  traces  over  the  same  time  windows  as  the  cross-correlation 
estimates. 
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The  suffix  A  on  filters  2  and  3  designates  a  modifi¬ 
cation  to  reject  rectilinear  motion  due  to  shear  phases  moving  in 
the  same  general  direction  as  the  compressional  phases. 

REMODE  5  is  an  experimental  polarization  filter  of 
greater  selectivity  than  REMODE  2  and  3.  It  is  designed  to 
attenuate  noise  having  nearly  rectilinear  polarization,  such  as 
signal-generated  noise,  which  is  not  effectively  rejected  by 
REMODE  2  or  3. 

As  a  test,  the  filters  were  applied  to  short-period 
LRSM  recordings  of  eight  large  teleseismic  events,  for  which 
reliable  focal  depth  estimates  were  available.  The  phases  pP 
or  sP  could  be  identified  with  confidence  from  three  of  the  events. 

Phases  which  were  clearly  visible  on  the  unfiltered 
seismograms  were  more  clearly  visible  on  the  normalized  filtered 
seismograms,  and  outstanding  on  the  un-normalized  filtered  seismo¬ 
grams.  All  of  the  depth  phases  which  were  identifiable  on  the 
filtered  seismograms  were  alio  identifiable  on  the  unfiltered 
seismograms. 

An  example  of  the  application  of  various  REMODE 
filters  is  shown  in  Figure  20. 

None  of  the  filters  effectively  enhanced  weak  signals. 
A  weak  sP  phase,  on  the  unfiltered  seismograms  of  the  Peru  event, 
was  slightly  more  evident  on  the  filtered  records.  A  tentative 
sP  phase  for  a  Bonin  Islands  event,  visible  on  several  of  the 
unfiltered  records  as  a  wavelet  broader  than  the  REMODE  window, 
was  degraded  by  un-normalized  filtering  and  unimproved  by  normal¬ 
ized  filtering. 
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If  no  trace  of  a  depth  phase  could  be  recognized 
across  the  suite  of  unfiltered  seismograms,  then  it  was  not 
recognizable  on  the  REMODE-processed  seismograms.  No  sP  phase 
could  be  found  on  any  seismograms  for  a  Colombia  event,  nor 
could  any  of  the  depth  phases  be  found  for  three  other  events. 

The  PcP  phase  could  not  be  identified  for  several  events.  Al¬ 
though  such  phases  could  be  presumed  to  exist,  there  is  no  reason 
to  assume  that  they  contain  enough  energy  to  be  detectable. 

2.  Feasibility  of  Linear  Polarization  Measurements  for 

Detecting  and  Measuring  Seismic  Body  Waves 

The  objective  of  this  study  was  to  compare  at  tele- 
seismic  distances  the  detectability  of:  P  on  the  vertical  component 
with  polarization  measurements. 

The  vertical  and  radial  component  of  a  strong  tele- 
seismic  event  was  added  to  noise  so  that  the  S/N  ratio  was  known 
a  priori.  Vertical  component  measurements  and  polarized  vertical 
component  measurements  were  compared  on  twenty  different  noise 
samples  at  approximately  twenty  different  noise  levels.  In 
addition,  a  theoretical  evaluation  was  made  assuming  Gaussian 
linearly  additive  noise.  The  sample  mean  amplitude  was  based  on 
twenty  samples.  On  the  assumption  that  the  noise  and  signal  were 
randomly  phased,  the  estimate  of  the  signal  amplitude  was  the 
square  root  of  the  difference  of  the  squared  signal  plus  noise 
and  squared  noise.  The  signal-to-noise  ratio  was  computed  Ly 
dividing  the  above  sampled  mean  of  signal  by  that  of  noise. 
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Results  of  the  evaluation  for  zero  lags  are  shown 
on  Figure  21  .  For  signals  with  input  S/N  greater  than  18  db, 
no  bias  in  the  amplitude  measurement  occurred  and  an  average 
gain  in  S/N  ratio  of  9  db  was  obtained  without  signal  distortion. 

For  input  S/N  greater  than  18  db,  a  second  pass  of  the  operator 
resulted  in  further  enhancement  of  several  db.  For  input  S/N 
from  9  to  18  db,  gains  were  obtained  from  0  db  to  9  db,  respect¬ 
ively;  but  at  the  expense  of  large  bias  in  amplitude  measurements 
as  shown  at  the  bottom  of  Figure  n  .  For  inputs  less  than  a  9  db 
threshold,  the  polarization  measurements  showed  a  loss  in  S/N 
ratio. 

Similar  results  were  obtained  using  the  even  part 
of  the  correlation  coefficient  up  to  +  second  of  lag.  Using 
lags,  input  greater  than  18  db  was  required  for  unbiased  ampli¬ 
tude  measurements,  but  a  gain  of  14  db  resulted.  With  lags, 
the  detection  threshold  was  near  5  db. 

The  gain  obtained  in  measuring  the  linearly  polarized 
component  on  Z  and  R  was  analyzed  theoretically,  assuming  the  signal 
is  linearly  polarized  and  the  ambient  noise  was  Gaussian  and  linearly 
additive.  In  order  to  evaluate  the  expected  gain  of  the  linearly 
polarized  measurement  compared  to  the  vertical  component  measurement, 
the  signal-to-noise  ratio  of  cross  correlated  energy  fluctuations 
of  the  two  components  was  compared  to  the  signal-to-noise  ratio  of 
the  energy  fluctuations  of  the  vertical  component.  It  was  found 
that  this  measure  of  gain  depends  only  on  the  noise  power,  noise 
coherency,  and  ratio  cf  horizontal  to  vertical  component  amplitude 
of  the  signal. 
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Figure  21.  Polarization  Measurements 
Using  REMODE  Filters 
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Based  on  the  test  for  polarization  measurements  used, 
no  significant  detection  gain  can  be  expected  uping  polarization 
criteria  to  detect  weak  teleseismir  P-vaves;  they  can  be  used  to 
isolate  and  measure  stronger  events.  For  Gaussian  noise,  thj 
measurements  of  linear  polarization  should  be  useful  at  epicentral 
distances  lees  than  40  for  which  several  db  gain  can  be  expected. 

3.  REMODE  Signal-to-Noi je  Tests  in  Polarized  Noise 

Much  of  the  noise  following  the  P  phase  had  frequency 
content  and  polarization  so  similar  to  signal  that  REMODE  proces¬ 
sing  did  not  isolate  weak  signals  effectively. 

The  REMODE  5  processor  xb  designed  to  pass  input 
motion  which  is  rectilinear  to  a  specified  degree  in  a  specified 
range  of  directions.  Its  sensitivity  to  rectilinear  motion  can 
be  controlled  by  variation  of  parameters. 

Short-period  test  tapes  were  generated  to  compare 
different  REMODE  5  processes  at  low  S/N  ratios.  The  test  tapes 
consisted  of  a  P-phase  embedded  in  its  own  signal-generated  noise, 
to  simulate  a  depth  phase  in  its  typical  noise  environment. 

A  tele&eism  from  the  Sea  of  Okhotsk  was  selected 
which  had  a  strong,  well-recorded  I’-phase  at  one  cf  the  LRSM 
stations.  The  first  four-cycle  segment  of  the  P  phase  on  the  Z 
component  was  multiplied  by  a  factor  and  added  to  both  the  Z  and 
R  components  of  the  noise  sample.  The  same  data  were  used  to 
represent  signal  in  shear  noise,  by  reversing  the  polarity  of 
the  radial-component  noise  before  addition  of  signal. 
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Figure  24. 
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6,500 

2,20 

1.5 

.5,1 

-10 

to 

-3  dB 

34 

13,500 

2,20 

.75 

.5,1 

-10 

to 

-3  dB 

22 

- 

2,20 

3. 

.5,1 

0 

to 

+  6  dB 

40 

- 

2,20 

1.5 

2. ,  1 

-10 

to 

-3  dB 

15 

- 

(REMODE 

1.5 

.5,1 

-10 

to 

-3  dB 

2.9 

230 

3A 

Table  5.  REMODE  Comparisons  for  Detection  Thresholds 
and  S/N  Improvements 


Figure  22  shows  radial-component  output  S/to  ratio 

estimates  p  as  functions  of  "true"  input  S/to  ratios  R,  for 

several  different  processors.  The  "input"  curve  represents 

input  S/to  estimated  in  the  same  way  as  output  S/to,  to  show  the 

bias  in  the  S/to  estimator  P.  Ideally  the  "input"  curve  would 

be  a  straight  line  of  slope  1.0.  Actually  its  slope  increases 

with  decreasing  S/to  because  of  cancellation  between  signal  and 

noise  of  opposite  polarity.  Figure  23  shows  output  S/to  estimates 

normalized  to  the  input  S/to  estimate,  and  represents  the  S/to 

improvement  or  gain  Y  of  each  processor  as  a  function  of  the 

jK 

"trite"  input  S/to  ratio  R.  The  processors  having  large  M  severely 
attenuate  shear  noise,  providing  large  S/N  improvement  for  input 
S/to  ratios  greater  than  1.0.  At  low  input  S/to  ratios,  trans- 
verse-polarized  noise  overrides  signal  polarization. 

Figures  24  and  25  show  the  output  S/N  ratios  and  S/N 
improvement  ratios  in  compressional  noise  for  the  radial  component. 
Output  S/N  ratios  for  most  of  the  processors  are  orders  of  magni¬ 
tude  smaller  than  they  were  for  shear  noise,  because  compressional 
noise  is  not  strongly  attenuated  by  polarization  filtering.  REMODE 
gain  (Figure  25)  increases  noticeably  with  increasing  input  S/N 
for  the  more  selective  processors.  The  negative  slopes  at  low 
S/to  ratios  are  due  to  differences  in  bias  between  output  S/to  es¬ 
timates  and  the  input  S/to  estimates  to  which  they  are  normalized. 

The  advantages  of  changes  in  IJMODE  selectivity  were 
most  evident  at  high  input  S/tl  ratios.  Although  the  processors 
differed  markedly  in  their  ability  to  improve  strong  signal,  they 
differed  much  less  in  their  ability  to  detect  signal  at  low  S/to 
ratios.  This  is  shown  in  Table  5,  which  summarizes  the  detection 
thresholds  and  maximum  S/N  improvement*  for  several  different 
processors. 
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Conclusions  regarding  the  use  of  REMODE  5  in  the 
detection  and  identification  of  depth  phases  include  the 
following. 

Increase  in  REMODE  selectivity  are  more 
effective  for  isolating  visible  signals  than  for  revealing 
obscure  signals. 

2.  There  exists  an  input  S/N  level  below  which 
REMODE  operators  will  provide  little  or  no  signal  enhancement. 

limit  is  below  the  S/N  level  at  which  signal  is  recognizable 
on  the  unprocessed  records. 

3.  Signal  enhancement  in  the  presence  of  compres¬ 
sions!  noise  does  not  necessarily  improve  signal  detectability  at 
low  S/N  levels,  because  nearby  compressional  noise  phases  larger 
than  the  signal  can  be  enhanced  even  more  than  the  signal. 


G.  Energy  Fluctuationa  in  Seismic  Noise 

In  recent  years,  statistical  methods  have  been  success¬ 
fully  used  to  predict  the  change  in  ambient  seismic  noise  power 
sensed  by  burying  a  seismometer  beneath  the  surface  of  the  earth. 
Observations  of  seismic  noise  has  indicated  a  large  standing  wave 
or  isotropic  component  composed  of  the  admixture  of  many  propaga¬ 
tion  modes.  Statistical  theories  such  as  equipartition  of  energy 
have  been  reasonably  successful  for  deriving  the  excitation  of 
the  propagation  modes. 

Since  such  models  are  classically  described  by  waves 
from  a  zero-mean  Gaussian  population,  the  purpose  of  this  study 
was  to  tentatively  assume  that  the  seismic  noise  is  Gaussian  and 
to  try  to  reject  the  hypothesis  by  measuring  the  relative  fre¬ 
quency  of  occurrence  of  specific  particle  energies  sequentially 
observed  on  samples  of  seismic  noise.  As  control  for  the  study, 
the  same  frequency  of  occurrence  observations  were  made  of  a 
thermal  noise  sample  recorded  on  tape  from  a  Gaussian  noise 
generator.  The  energy  envelope  of  Gaussian  nol*e  is  theoret¬ 
ically  described  by  a  Boltzman  or  exponential  probability  function. 
Thus  the  fit  of  a  straight  line  to  the  log  of  the  relative  fre¬ 
quency  of  occurrence  is  a  test  of  the  Gaussian  hypothesis.  Two 
tests  were  used;  one  was  the  chi-square  test  of  the  deviation 
from  the  least  squares  straight  line;  the  other,  more  powerful 
test  compared  the  variance  of  the  seismic  noise  and  that  of  the 
Gaussian  noise  generator. 

Three  hours  of  seismic  noise  was  selected  which  appeared 
to  be  typical  of  the  normal  ambient  background.  The  noise  was 
measured  at  a  depth  of  7,452  feet  in  a  deepwell  near  Apache,  Okla¬ 
homa.  The  near  surface  layering  is  described  as  high  velocity 
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limestones  of  approximately  6  km/sec  overlying  an  igneous  basement 
complex  of  velocities  in  the  neighborhood  of  5  km/sec. 

As  shown  on  Figure  26  for  seismic  noise,  the  relative 
frequency  of  occurrence  is  given  by  the  vertical  axis,  and  the 
class  interval  in  order  of  increasing  energy  is  shown  on  the 
horizontal  axis.  On  each  plot,  the  center  frequency  and  effect¬ 
ive  bandwidth  of  the  filter  is  labeled.  On  Figure  27  ,  the  same 

information  is  shown  for  the  Gaussian  noise  generator.  Comparing 
Figure  26  with  Figure  27  ,  the  scatter  from  a  linear  trend  is  the 
same  overall  for  both  sets  of  data,  seismic  and  Gaussian  generator. 

Based  on  the  chi-square  test,  it  is  more  than  99%  prob¬ 
able  that  samples  at  all  of  the  frequencies  are  from  a  zero  mean 
Gaussian  population.  The  hypothesis  that  the  class  deviation  of 
the  seismic  noise  sample  and  the  thermal  noise  sample  are  from 
the  same  random  population  is  rejected  only  for  the  sample  at 
1.4  cps,  suggesting  that  a  non-Gaussian  component  may  be  included 
with  the  noise  in  this  band.  If  we  accept  the  Gaussian  hypothesis 
theoretical  interest  in  the  power  spectral  density  of  the  noise 
may  bear  ultimately  on  dissipation  mechanisms  and  may  possibly 
lead  to  analysis  of  geological  structure,  with  the  real  changes 
in  the  power  spectrum  in  a  region  depending  primarily  on  changes 
in  structure.  Of  more  practical  interest  is  that  modern  litera¬ 
ture  on  detection  and  filtering  is  most  meaningful  in  the  context 
of  seismic  signals  added  to  Gaussian  noise. 


Figure  26.  Relative  Frequency  of  Occurrence  (Ordinate)  Versus 

Energy  Level  (Abscissa) 


H.  The  Coherency  Analysis  of  Seismic  Noise 

In  the  coherency  analysis  of  seismic  noise  at  arrays, 
we  would  like  to  know  the  fraction  of  noise  which  is  correlated 
and  uncorre.lated  at  each  frequency,  the  linear  relations  between 
seismometer  outputs,  the  stationarity  of  the  noise,  and  other 
properties  in  an  effort  to  predict  how  well  optimum  filters  could 
be  made  to  work  if  computer  limitations  were  no  problem.  The 
theory  is  used  to  estimate  frequency  response  functions,  ordinary 
coherence  functions,  partial  coherence  functions,  and  multiple 
coherence  functions. 

The  parameters  being  estimated  may  be  qualitative  j.y 
interpreted  as  follows.  The  valve  of  a  frequency  response  func¬ 
tion  at  frequency  point  f  can  be  expressed  in  terras  of  gain  and 
phase.  The  gain  describes  the  change  in  amplitude  and  the  phase 
describes  the  phase  shift,  or  equivalently  the  time  delay,  ex¬ 
perienced  by  a  sinusoid  of  frequency  f  passing  through  the  linear 
time  invariant  system  characterized  by  the  frequency  response 
function.  The  ordinary  coherence  function  provides  a  measure 
of  the  degree  of  linear  relationship  between  two  variables  when 
they  are  considered  as  the  input  and  output  of  the  time  invariant 
linear  system.  An  alternative  equivalent  interpretation  is  that 
of  the  percentage  of  power  or  variance  which  is  common  to  the  two 
random  variables. 


The  partial  coherence  function  is  a  quantitative 
measure  of  linear  relationship  between  two  variables  when  they 
are  considered  as  being  related  by  a  linear  time  invariant  system 
and  linear  least  squares  estimates  of  other  measured  variables 
have  been  subtracted  out.  Under  appropriate  conditions,  the 


partial  coherence  between  two  variables  will  provide  a  more 
pertinent  measure  of  their  relationship  than  will  ordinary  coher¬ 
ence  functions.  The  multiple  coherence  function  provides  a 
measure  of  the  degree  of  linear  relationship  between  a  single 
record  considered  as  an  output  variable  and  several  other 
variables  considered  input  variables.  An  alternative  interpre¬ 
tation  of  multiple  coherence  is  the  percentage  of  power  or 
variance  which  is  accounted  for  via  linear  relationships  with 
all  the  input  variables  simultaneously. 

The  multiple  coherence  function  cm  assist  in  deter¬ 
mining  how  many  seismometer  outputs  in  an  array  should  be  proces¬ 
sed  to  properly  determine  the  seismic  .>oise  field.  The  partial 
coherence  functions  and  ordinary  coherence  functions  can  assist 
in  the  determination  of  where  pertinent  strong  linear  relation¬ 
ships  exist  from  seismometer  to  seismometer  to  help  direct  efforts 
eventually  intended  toward  noise  suppression  with  a  minimum  of 
computational  effort. 

Digital  computer  programs  were  de  reloped  and  sets  of 
numerical  experiments  were  made.  These  numerical  experiments 
illustrated  that  the  multiple  coherence  function  acts  basically 
as  predicted.  It  attains  a  low  value  when  there  are  less  seis¬ 
mometer  outputs  being  measured  than  there  are  underlying  noise 
traces,  and  it  increases  and  approaches  unity  as  the  number  of 
seismometers  increases  past  the  number  of  underlying  independent 
noise  traces. 

The  partial  coherence  function  is,  in  general,  more 
difficult  to  interpret  in  that  residual  coherences  can  exist 
between  linear  combinations  of  extraneous  noise.  However,  it 
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would  appear  that  if  some  assumptions  are  made  concerning 
the  character  of  the  incoherent  extraneous  noise  (e.g., 
constant  power  spectra  throughout  the  seismometer  array) 
then  one  should  be  able  to  predict  the  residual  partial 
coherence.  This  prediction  can  be  subtracted  out  such  that 

the  remainder  is  close  L'  zero. 

The  gain  and  phase  information  relating  the  seis¬ 
mometer  output  has  yet  to  be  exploited  in  detail,  but  it 
should  be  useful  in  determining  direction  and  speed  of  the 
noise  components.  Other  information  available  from  the 
program  is  the  conditional  power  spectral  density  functions. 

Thi*  x a  the  power  remaining  alter  linear  combinations  of  other 
variables  have  been  subtracted  out.  This  gives  additional 
data  to  provide  more  insight  into  the  structure  of  the  noise. 

As  an  example  of  the  kind  of  approach  that  might  be 
applied  to  an  array  using  coherency  theory,  the  Cumberland 
Plateau  seismic  array  given  in  Figure  28  was  examined.  The  sub- 
array  analyzed  contained  nine  seismometers  with  the  basic  data 
being  a  noise  sample  of  the  nine  traces,  40C  seconds  long,  taken 
at  5  points  per  second  with  50  lags.  A  first  consideration  wa 
the  determination  of  characteristics  of  this  noise  field  which 
related  the  amount  of  noise  power  in  seismometer  10  that  could 
be  predicted  or  accounted  for  by  taking  various  combinations 
of  other  elements  in  the  array.  Thus,  the  center  was  considered 
an  "output"  and  various  combinations  of  other  seismometers  were 
sequentially  added  au  "inputs"  to  determine  the  change  in  the 
multiple  coherence,  i  e.  \e  additional  proportion  of  power 
accounted  for  by  each  i  i. 
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Figure  29  gives  the  multiple  coherence  or  the  propor¬ 
tion  of  the  nci  se  power  in  10  which  was  accounted  for  by  3,  4, 
5,...,  9  component  seismometer  arrays,  respectively,  at  five 
different  frequencies.  The  outer  seismometers  were  added  in 
first  in  this  case.  We  note  that  much  of  the  power  (94.9%) 
in  the  .2  cps  band  is  already  accounted  for  by  the  first  three 
seismometers.  Figure  29  also  shows  the  rate  of  increase  in  the 
proportion  of  power  in  10  accounted  for,  as  the  number  of  seis¬ 
mometers  increases.  At  low  frequencies,  most  of  the  power  is 
accounted  for  by  the  time  seven  seismometers  are  present  while 
at  1.8  cps  the  full  set  accounts  for  only  60  percent  of  the 
noise  power.  By  examining  the  ordinary  coherence  as  given  in 
Table  6 ,  one  notices  that  seismometer  10  is  highly  coherent 
with  most  of  the  other  seismometers  tending  to  verify  that  10 
is  a  reasonable  candidate  for  the  "output"  variable  in  one  par¬ 
ticular  linear  system  represented  by  the  array.  However,  the 
high  ordinary  coherence  between  10  and  2,  8  or  6  at  low  frequen¬ 
cies  suggests  that  adding  t'nose  closer  seismometers  fjirst  would 
have  accounted  for  more  of  the  coherent  no.’ se  power.  In  fact 
it  can  be  seen  from  the  ordinary  coherence  in  Table  6  that 
adding  seismometer  8  into  the  array  first  would  have  produced 
a  two-element  array  which  accounted  for  between  30  and  98%  of 
the  noise  power,  an  arrangement  clearly  superior  to  the  first 
three-element  array  in  Figure  28.  Thus,  by  choosing  certain 
seismometers  in  the  array  one  could  improve  the  slope  of  the 
"proportion  of  noise  power  accounted  for"  curve  in  Figure  29. 
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In  Figure  30/  the  same  sample  (A)  is  considered  with 
the  difference  being  the  addition  of  the  close-in  seismometers 
first.  This  time  the  proportion  of  power  accounted  for  increases 
much  more  sharply  as  the  number  of  seismometers  increases  and  by 
the  time  six  are  present  the  noise  accounted  for  compares  favor¬ 
ably  with  the  previous  situation  when  all  nine  are  present.  In 
fact  more  than  90%  of  the  power  in  the  low  three  bands  is  accounted 
for  by  a  three-element  array,  consisting  of  the  three  close-in 
seismometers. 

The  data  in  both  Samples  A  and  B  indicate  that  the 
outer  elements  are  coherent  with  the  center  seismometer  10 
(Table  6) .  This  suggests  that  much  of  the  noise  power  in  the 
first  two  bands  could  be  eliminated  by  replacing  each  seismometer 
by  its  residual  trace  (i.e.,  with  the  effect  of  10  eliminated). 

With  this  in  mind  the  simple  linear  filter  relations  between  10 
and  the  other  seismometers  (given  in  Table  7)  wire  examined  and 
found  to  be  quite  consistent  over  the  two  racording  samples. 

This  shows  that  for  the  inner  seismometers  (See  Table  6)  between 
92  and  100%  of  the  noise  power  can  be  accounted  for  at  each 
seismometer  for  .2  cps  noise  and  70  through  90%  of  the  noise 
power  at  .6  cps  merely  by  subtracting  out  the  effect  of  10  and 
furthermore,  that  it  is  done  in  both  situations  with  the  same 
set  of  filter  coefficients.  This  still  leaves  plenty  of  degrees 
of  freedom  left  for  estimating  the  signal  and  gets  rid  of  the 
noise  with  the  same  set  of  filters  in  both  situations. 
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Figure  30.  Sample  A  —  Proportion  of  Noise  Power  Accounted  For 
With  Inner  Seiemometers  Added  First 


Sample  A  ( . 2  cps ) 


Real 

1.35 

1.16 

.94 

1.06 

Imaginary 

.56 

.20 

39 

.14 

Sample  B  (.2  cps) 


Real 

1.33 

1.19 

.91 

1  08 

Imaginary 

.54 

.19 

.38 

.13 

Sample  A  (.6  cps) 

Real 

1.24 

.96 

.89 

.66 

Imaginary 

.09 

-.14 

.12 

-.60 

Sample  B  (.6  cps) 

Real 

1.16 

.98 

.81 

.63 

Imaginary 

.11 

-.11 

.13 

-.61 

Table  7.  Linear  Filter  Relation  Between  Sub-Array  Element 

and  Center  Element  (10) 


I.  Analysis  of  Seismic  Noise  at  the  Yellowknife  Array, 

Canada 

From  data  supplied  by  the  United  Kingdom  Atomic  Energy 
Authority  (UKAEA)  a  quiet  day  noise  analysis  wsa  performed  for 
the  Yellowknife  array  (YKC) ,  Canada.  The  data  was  supplied  in 
the  form  of  a  twenty-fdur  channel  magnetic  tape  recording  of  the 
amplified  outputs  of  short-period,  vertical-component,  Willmore 
Mark  II,  seismometers.  Beth  an  original  recording  and  a  trans¬ 
cription  were  supplied.  The  UKAEA  transcription  process  was 
analyzed  for  fidelity  by  examining  several  noise  samples,  taken 
at  vhe  same  time  intervals,  from  the  original  and  transcription 
tapes.  Although  the  noise  samples  of  the  transcription  looked 
very  much  like  the  original  recordings,  power  spectral  analysis 
indicated  generally  high  fidelity  only  below  about  1.5  cps. 

Above  1.5  cps  system  generated  noise  led  to  poorer  reproduction. 

Analysis  of  the  noise  field  is  advantageously  per¬ 
formed  in  terms  of  the  power  in  the  field  as  a  function  of  fre¬ 
quency  and  vector  wave-number  i.e.,  3-D  powdr  spectral  analysis 
where,  for  any  frequency,  the  power  is  given  as  a  function  of 
tho  vector  wave-number,  k.  However,  the  relative  spacing  of  the 
seismometers  at  YKC  caused  severe  aliasing  leading  to  uninter- 
protable  spectra.  Thus  an  attempt  was  made  to  arrive  at  a 
meaningful  picture  of  the  vertical  component  of  the  noise  field 
by  calculating  auto-  and  cross-correlation  functions,  single 
channel  power  spectra,  crcss-power  spectra,  end  coherencies. 

The  array  geometry  and  total  system  response  are  shown  in  Figures 
31  and  32,  respectively. 
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The  field  appears  tc  be  divided  into  coherent  and  in¬ 
coherent  (random)  noise  where  the  coherency  depends  on  frequency, 
seismometer  spacing  and  direction.  The  noise  below  about  .4  cps 
generally  exhibits  the  highest  coherency  with  a  peak  of  90%  for 
a  2.5  km  spacing  along  the  B  arm  of  the  array  and  decreases  with 
increasing  frequency  and  seismometer  separation.  The  peak  co¬ 
herency  for  it  2.5  km  separation  along  the  R  seismometers  is  72%, 
and  drops  more  rapidly  with  increasing  spacing  than  along  the  B 
arm.  Instrumental  nois*  appears  to  contribute  to  incoherent  power 
at  all  frequencies  and  to  coherent  power  at  several  frequencies. 

In  particular,  instrumental  noise  was  found  to  contribute  such 
large  amounts  of  incoherent  as  well  as  coherent  power  as  to  make 
analysis  above  about  1.5  cps  impossible. 

The  peak  microseismic  energy  was  found  to  occur  at  a 
period  of  4.5  seconds,  at  abrjt  20  db  above  the  1  cps  noise  power 
as  seen  on  spectra  which  were  uncorrected  for  the  system  response, 
'’’he  4.5  second  peak  appears  to  consist  of  a  wave  or  group  of  waves 
traveling  at  about  3.6  km/sec.  in  an  easterly  direction  with  a 
small  northerly  component.  This  dominant  micro seismic  energy  hits 
the  B  seismometers  almost  broadside  thus  accounting  for  the  apparent 
directional  dependence  of  the  low  frequency  coherence.  Corrected 
for  system  response,  the  peak  occurs  at  about  8.0  seconds  at  70-80 
db  above  the  1  cps  power. 

It  is  obvious  from  the  above  discussion  that  the  dynamic 
range  of  the  instrumentation  must  be  known  in  deta'l  and  that 
analysis  of  the  noise  field  above  about  1.5  cps  requires  either 
increasing  the  overall  dynamic  range  or  eliminating  the  dominant 
low  frequency  energy  before  sensing  and  recording.  The  geometry 
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of  the  array  determines  the  spatial  spectral  window  and  thua  the 
power  in  the  noise  field  that  is  "seen”  in  k-space.  That  ia, 
analogous  to  the  Nyquist  folding  frequency,  the  spatial  sampling 
determines  the  spatial  folding  vector  wave-number  and  also  the 
shape  and  position  (in  k-space)  of  the  sidelobea.  For  YKC  the 
seismometers  were  too  far  apart  causing  severe  alissing  thua 
making  a  k-space  noise  analysis  impossible.  The  regular  spacing 
of  the  array  is  desirable  since  this  causes  shape  peaks  in  the 
side  lobes,  but  a  closer  spacing  (at  most  .5  km  apart)  is  neces¬ 
sary  in  order  to  increase  the  spatial  folding  wave-number  to  the 
point  where  a  meaningful  k-apace  analysis  is  possible. 

Figures  33  and  34  are  illustrative  of  the  type  of  analysis 
made  for  this  study.  Figurr  33  shows  the  original  and  associated 
transcription  of  one  sample  taken  at  the  Rl  seismometer.  Figure 
34  shows  low  resolution  auto-spectra  computed  for  the  Rl  seis¬ 
mometer  for  the  sample  given  in  Figure  33. 
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Figure  33.  Yellowknife  Array  Noiee  Sample* 
Original  and  Transcription 


*  . 


iJ 


.  finite  Fourier  Transform  Theory  and  It*  Application  to 


the  Computation  of  Convolution*,  Correlations,  and  Spectra 


Practical  end  computational  aspect*  of  the  theory  of 


Fourier  Transforms  have  been  examined  in  connection  with  various 


SDL  research  analyses.  These  efforts  have  resulted  in  a  set  of 
programs  for  performing  operations  on  time  series  based  on  the 
Cooley-Tukey  hyper-rapid  Fourier  transform  method.  Using  this 
method,  computations  on  seismic  array  data  such  as  the  calculation 
of  convolutions,  correlations,  spectra,  and  digital  filters  have 
been  speeded  up  by  factors  of  three  or  four  and  sometimes  even  ten. 
Following  is  a  brief  description  of  the  procedures  followed  and 
the  results  obtained  from  these  procedures. 


1.  Finite  and  Discrete  Fourier  Transforms 


In  the  case  of  continuous  data  of  infinite  length,  the 
Fourier  transform  pair  is  the  direct  transform  and  the  inverse 


transform.  Sometimes  the  direct  transform  is  written  with  a  factor 


of  1  in  front  on  the  integral  and  the  inverse  with  a  factor  of  1/2  tt 
Quantities  of  interest,  such  as  spectra,  etc.,  involve  magnitudes 
or  squares  of  on*  transform  and  the  factor  must  be  inserted  or 
taken  out,  depending  on  which  definition  is  used,  to  preserve  the 
proper  dimensions. 

Two  drawbacks  of  these  definitions  for  digital  com¬ 
putations  are  apparent:  First,  the  integrals  must  be  approximated 
by  sums  in  the  digital  computer,  which  implies  that  both  trans¬ 
forms  Involve  sampled  variables.  Second,  the  infinite  limits  on 
the  sums  are  impossible  in  practice.  Clearly,  these  sums  must 
be  truncated,  as  they  do  not  in  general  converge  over  a  finite 
interval.  As  a  result,  Fourier  transforms  as  such  are  never  really 
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computed  by  a  digital  computer.  Ir  stead,  the  complex  samples 
of  a  direct  transform  are  approximated  by  the  coaine  and  aine 
coefficients  of  Fourier  aeriea  repreaentation  of  the  input  data. 

If  N  aamplea  of  the  data  are  taken  at  equally 
apaced  intervala  At  ■  T/N,  the  integrala  become  sums,  and  the 
frequency  aum  goea  from  DC  to  the  folding  frequency.  It  waa 
found  that  a  great  deal  of  aymmetry  between  the  two  tranaforma 
could  be  preaerved  if  the  aum  ia  summed  up  to  N-l.  Redundant 
pointa  in  the  apectrum  are  included  (aince  the  tranaforma  are 
periodic)  but  the  computational  procedurea  are  aimplified. 

It  waa  ahown  that  the  aet  of  direct  Fourier  trans¬ 
form  pointa,  between  DC  and  the  folding  frequency,  contained  the 
same  amount  of  information  as  the  renl  data  series,  which  sug¬ 
gested  that  the  existence  of  one  transform  should  imply  the 
existence  of  the  other. 

2.  Two-and-Three-Dimensional  Fourier  Transforms 

Two-  and  three-dimensional  direct  Fourier  tranaforma 


are  seen  to  be 


A(k1#k2)  - 


Vi  V1 


N1N2  j  -0 


j2-o 


x(J1»j2)  w1”Jlkl  w2“^2k2 


A(klfk2,k3) 


Nx-1  N2-l  N3-l 


l  l 


N1N2N3  jx-0  j2-0  j3-0 


x(jrj2>j3)  w  “31K1. 


.w2":)2k2  w3”^3k3 


.  - 


* 

S' 

\  ' 


A  ‘ 


By  separating  and  breaking  up  these  equations,  it  was 
calculated  that  +  N2  one-dimensional  transforms  are  required  to 
compute  the  single  two-d tmensional  transform,  and  that  N,N2  one¬ 
dimensional  transforms  and  two-dimensional  transforms  are  needed 
to  compute  the  single  three-dimensional  transform. 

3.  Speed-Up  of  Transform  Computing  Time 

It  was  shown  that  the  process  of  transforming  is 
equivalent  to  matrix  multiplication  by  a  matrix  W,  which  preserves 
"length"  between  two  domains.  The  Cooley-Tukey  **  Jthod  factors  the 
W  matrix,  if  its  order  is  a  power  of  two,  into  L  +  1  sparse  matrices, 
where  L  is  the  power  of  two.  Multiplying  L  +  1  tiroes  by  these  sparse 
matrices  can  in  some  cases  reduce  the  computing  time  by  many  tens  of 
times. 

4.  High-Speed  Correlations  and  Convolutions 

By  computing  Fourier  transforms  with  a  finite  Fourier 
series-like  method  an  important  condition  is  put  on  the  time  series. 
As  in  regular  Fourier  series  the  input  is  assumed  to  be  periodic 
with  period  T  and  the  integrals  or  sums  are  computed  over  a  single 
period.  There  is  also  the  effect  of  cutting  off  the  spectrum  at 
the  folding  frequency.  Sines  and  cosines  of  finite  wavelength  will 
repeat  again  outside  the  region  of  interest.  This  fzct  in  itself 
is  not  bothersome  but  becomes  a  serious  complication  in  the  computa¬ 
tion  of  convolutions  and  correlations.  Convolutions  and  correla¬ 
tions  as  usually  computed  assume  the  time  series  to  be  zero  outside 
the  region  of  interest.  Therefore,  the  integrals  or  sums  in  com¬ 
puting  them  are  summed  out  only  o/er  the  non-zero  terms.  When 
multiplying  together  two  finite  Fourier  transforms  (or  the  complex 
conjugate  of  one  times  the  other)  the  periodicity  of  the  time 


series  mean*  that  elements  which  have  been  shifted  pa  it  the  end 
of  a  period  reappear  at  the  beginning.  This  process  is  called 
circular  convolution  or  correlation  and  its  effects  are  un¬ 
avoidable  when  straightforwardly  computing  lagged  products  with 
finite  Fourier  transforms. 

Circular  convolution  is  written: 

R^(t)  -  Tjj1  x.(t  )  x.(t+  t)  (3) 

where  x  (t  +  T)  ■  x  (t)  for  all  m. 
m  m 

It  was  shown  that  this  kind  of  correlation  is  equal 
to  the  transform  of  the  absolute  product  of  the  two  finite  trans¬ 
forms. 

On  the  other  hand  the  transient  correlation  for 
positive  lags  is  defined  by  the  following: 

Rii(t)  “  T"£_t  xi(T)  xl(t  +T  }  (4) 

J  f  «0  ^ 

where  the  upper  limit  on  the  sum  simulates  tho  desired  zeros  in 

the  time  series  outside  the  region  of  interest.  The  finite 

T 

Fourier  transform  of  this  R  is  b«us  not  the  product  of  the  two 
individual  transforms.  However,  by  filling  zeros  into  the  second 
half  of  each  vlata  series  and  computing  their  transforms  out  to 
twice  their  actual  length,  a  good  estimate  of  the  spectrum  may  be 
obtained.  In  addition,  the  negative  lags  in  the  correlation  appear, 
thus  giving  a  more  mathematically  satisfying  result. 

Transient  correlations  for  100%  lags  were  s^own  to  be 
computed  by  forming  tho  absolute  product  of  two  transforms,  each 
computed  out  to  twice  the  length  of  the  original  data  series  with 
zeros  filled  into  the  second  halves. 
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Non-circular  or  transient  convolutions  were  also 
computed  in  much  the  same  way,  except  that  the  transforms  had  to 
be  computed  out  to  a  length  equal  to  the  sum  of  the  lengths  of 
the  time  series  and  the  filter,  with  the  appropriate  number  of 
zeros  filled  into  each.  The  convolution  theorem  was  proved  in 
the  same  fashion  and  were  computed  by  forming  the  product  of  the 
two  transforms,  each  computed  out  to  a  length  equal  to  their  sum 
with  zeros  filled  into  the  extra  lengths. 


K.  Perturbation  Theory  for  the  Inveraion  of  Travel 

Time  Data 

Estimation  of  the  internal  elastic  structure  of  the 
earth  from  body  wave  travel  times  ordinarily  requires  both  arrival 
time  data  and  the  slope  of  a  curve  fitted  to  the  observed  travel 
time.  There  arc  well-known  difficulties  in  fitting  such  curves, 
particularly  in  cases  where  the  observed  travel  times  suggest  a 
multivalue  travel  time  function.  We  have  considered  methods  which 
make  use  of  later  arrival  Information  in  order  to  account  for 
possible  multiplicities  in  the  travel  time  curves.  An  additional 
advantage  to  be  gained  from  multiple  phase  inversion  is  the  re¬ 
duction  in  the  number  of  station  observations  necessary  to  deter¬ 
mine  earth  structure.  The  usefulness  of  the  bo 'y  phase  amplitudes, 
in  addition  to  the  travel  time  information,  was  also  demonstrated, 
since  the  amplitude  information  in  the  region  of  multiplicities 
can  be  used  to  confirm  the  actual  presence  of  a  multiplicity,  and 
in  addition,  the  nature  of  the  velocity  variation  giving  rise  to 
the  multiplicity. 

It  is  clear  that  an  effective  means  of  including  later 
arrival  data,  as  well  as  data  from  a.*.y  number  of  different  phases, 
is  through  use  of  an  iteration  scheme.  Such  an  approach,  which 
does  not  require  travel  time  3lope  information,  car;  incorporate 
constraints  and  information  from  sources  other  than  the  travel 
time  data.  In  particular,  simultaneous  inclusion  of  amplitude 
and  travel  time  information  in  an  inversion  procoee  using  both 
perturbation  and  iteration  can  be  accomplished,  although  the 
amplitude  information  would  sewc  only  as  a  condition  of  accept¬ 
ance  for  structure  which  otherwise  satisfied  the  travel  time 
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data  in  the  least  squares  sense.  What  is  required  for  the 
implementation  of  an  iteration  procedur*  is  a  general,  sys¬ 
tematic  and  accurate  means  of  generating  the  correction  to  be 
applied  to  the  structure  at  each  iteration  stage.  A  properly 
conceived  variational  treatment  has  been  used  to  supply  this 
analytical  informaticn. 

The  inference  of  earth  structure  from  body  wave  d*ta 
har  been  developed,  using  procedures  similar  to  those  used  in 
determining  elastic  and  tnelastic  earth  structure  from  surface 
wave  amplitude  and  dispersion.  The  method  developed  iteratively 
perturbs  an  initial  velocity  distribution  (which  may  be  shown 
to  be  in  agreement  with  surface  wave  dispersion  data  over  the 
region  in  question)  until  the  theoretical  travel  time  function 
is  in  agreement,  in  the  least  squares  cense,  with  the  observed 
data.  The  method  does  not  require  any  estimate  of  the  curve 
slope  or  multiplicities,  nor  does  it  require  a  particularly 
dense  and  correlate  data  coverage  over  the  entire  distance  range. 
Several  phases  can  be  used  simultaneously. 
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